/*
 * case.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#include "case_2D3t_SGA.h"
// MY HEADERS
#include"mros-v42.h"

namespace case_2D3t_SGA
{
global glob;

// *******************************************************************************************  //
// *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY ***
// *******************************************************************************************  //

#define Sqrt sqrt
#define Power pow

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double n;
	double J;
	double t;
	double tPrim;
	double tPrim2;
	double U;
	double beta;

	unsigned int L; //8
};

struct variables {
	double barchiAB;
	double simchiAB;
	double chiS;
	double chiT;
	double chiSPrim;
	double chiTPrim;
	double DeltaS;
	double DeltaT;
	double simDeltaB;
	double mu;
	double mAF;
	double d;
    double Energy; //13
 };

/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

inline void clear(parameters &p) {
	p.n = p.J = p.t = p.tPrim = p.tPrim2 = p.U = p.beta = 0.; p.L = 0;
}

inline void clear(variables &v) {
	v.barchiAB = v.simchiAB = v.chiS = v.chiT = v.chiSPrim = v.chiTPrim = v.DeltaS = v.DeltaT = v.simDeltaB = v.mu = v.mAF = v.d = 0;
	v.Energy = 0;
}

inline void addTo(const variables &q, variables &v, double d=1.) {
	v.barchiAB += q.barchiAB*d;
	v.simchiAB += q.simchiAB*d;
	v.chiS += q.chiS*d;
	v.chiT += q.chiT*d;
	v.chiSPrim += q.chiSPrim*d;
	v.chiTPrim += q.chiTPrim*d;
	v.DeltaS += q.DeltaS*d;
	v.DeltaT += q.DeltaT*d;
	v.simDeltaB += q.simDeltaB*d;
	v.mu += q.mu*d;
	v.mAF += q.mAF*d;
	v.d += q.d*d;
	v.Energy += q.Energy*d;
}

inline void copyTo(variables &v, const variables &q, double d=1.) {
	v.barchiAB = q.barchiAB*d;
	v.simchiAB = q.simchiAB*d;
	v.chiS = q.chiS*d;
	v.chiT = q.chiT*d;
	v.chiSPrim = q.chiSPrim*d;
	v.chiTPrim = q.chiTPrim*d;
	v.DeltaS = q.DeltaS*d;
	v.DeltaT = q.DeltaT*d;
	v.simDeltaB = q.simDeltaB*d;
	v.mu = q.mu*d;
	v.mAF = q.mAF*d;
	v.d = q.d*d;
	v.Energy = q.Energy*d;
}

inline void copyTo(parameters &v, const parameters &q) {
	v.n = q.n;
	v.J = q.J;
	v.t = q.t;
	v.tPrim = q.tPrim;
	v.tPrim2 = q.tPrim2;
	v.U = q.U;
	v.beta = q.beta;
	v.L = q.L;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters  p.n  p.J  p.t  p.tPrim  p.U  p.beta  p.L
	variables v; // variables  v.chiAB  v.chiS  v.chiT  v.barDelta  v.simDelta  v.mu  v.mAF  v.d
	double n, d, m, mu; // contractions

	// Abbreviations - dependend on kx, ky
	double gammak;
	double etak;
	double gammaPrimk;
	double gamma2Primk;

public:
	double E[4];
	double dEdlbarchiAB[4];
	double dEdlsimchiAB[4];
	double dEdlchiS[4];
	double dEdlchiT[4];
	double dEdlchiSPrim[4];
	double dEdlchiTPrim[4];
	double dEdlDeltaS[4];
	double dEdlDeltaT[4];
	double dEdlsimDeltaB[4];
	double dEdln[4];
	double dEdlmAF[4];

	double gt, gs, gtPrimMin, gtPrimMax, gtPrim2Min, gtPrim2Max;
	 double dgtdm, dgtdd, dgtdn;
	 double dgsdm, dgsdd, dgsdn;
	 double dgtPrimMindm, dgtPrimMindd, dgtPrimMindn;
	 double dgtPrimMaxdm, dgtPrimMaxdd, dgtPrimMaxdn;
	 double dgtPrim2Mindm, dgtPrim2Mindd, dgtPrim2Mindn;
	 double dgtPrim2Maxdm, dgtPrim2Maxdd, dgtPrim2Maxdn;
	double lbarchiAB, lsimchiAB, lchiS, lchiT, lchiSPrim, lchiTPrim, lDeltaS, lDeltaT, lsimDeltaB, ln, lmAF;

	Func(const parameters &p, const variables &v) {
		set(p, v);
		//refresh();
	}

	void set(const parameters &p, const variables &v) {
		copyTo(this->p, p);
		copyTo(this->v, v);
		refresh();
	}

	void set(const parameters &p) {
		copyTo(this->p, p);
		refresh();
	}

	// Gutzwiller
	double fgt();
	double fgs();
	double fgBDelta();
	double fgm();
	double fgtPrimMin();
	double fgtPrimMax();
	double fgtPrim2Min();
	double fgtPrim2Max();

	double fdgtdm();
	double fdgtdd();
	double fdgtdn();

	double fdgsdm();
	double fdgsdd();
	double fdgsdn();

	double fdgtPrimMindm();
	double fdgtPrimMindd();
	double fdgtPrimMindn();
	double fdgtPrimMaxdm();
	double fdgtPrimMaxdd();
	double fdgtPrimMaxdn();

	double fdgtPrim2Mindm();
	double fdgtPrim2Mindd();
	double fdgtPrim2Mindn();
	double fdgtPrim2Maxdm();
	double fdgtPrim2Maxdd();
	double fdgtPrim2Maxdn();

    // Energy
    double E1(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double E3(double gammak, double gammaPrimk, double gamma2Primk, double etak);

    double dE1dl_barchiAB(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE1dl_simchiAB(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE1dl_chiS(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE1dl_chiT(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE1dl_chiSPrim(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE1dl_chiTPrim(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE1dl_n(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE1dl_mAF(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE1dl_DeltaS(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE1dl_DeltaT(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE1dl_simDeltaB(double gammak, double gammaPrimk, double gamma2Primk, double etak);

    double dE3dl_barchiAB(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE3dl_simchiAB(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE3dl_chiS(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE3dl_chiT(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE3dl_chiSPrim(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE3dl_chiTPrim(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE3dl_n(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE3dl_mAF(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE3dl_DeltaS(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE3dl_DeltaT(double gammak, double gammaPrimk, double gamma2Primk, double etak);
    double dE3dl_simDeltaB(double gammak, double gammaPrimk, double gamma2Primk, double etak);

	// DEBUG
	void printLambda() {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e (l)\n", lbarchiAB, lsimchiAB, lchiS, lchiT, lchiSPrim, lchiTPrim, lDeltaS, lDeltaT, lsimDeltaB, ln, lmAF);
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e (g)\n", gt, gs, gtPrimMin, gtPrimMax, gtPrim2Min, gtPrim2Max);
	}

	// REFRESH
	void refresh() {
		n = p.n;
		d = v.d;
		mu = v.mu;
		m = v.mAF;

		// obliczam funkcje Gutzwillara i poszczegolne pochodne
		gt = fgt();
		gs = fgs();
		gtPrimMin = fgtPrimMin();
		gtPrimMax = fgtPrimMax();
		gtPrim2Min = fgtPrim2Min();
		gtPrim2Max = fgtPrim2Max();

		dgtdn = fdgtdn();
		dgtdd = fdgtdd();
		dgtdm = fdgtdm();

		dgsdn = fdgsdn();
		dgsdm = fdgsdm();
		dgsdd = fdgsdd();

		dgtPrimMindm = fdgtPrimMindm();
		dgtPrimMindd = fdgtPrimMindd();
		dgtPrimMindn = fdgtPrimMindn();
		dgtPrimMaxdm = fdgtPrimMaxdm();
		dgtPrimMaxdd = fdgtPrimMaxdd();
		dgtPrimMaxdn = fdgtPrimMaxdn();

		dgtPrim2Mindm = fdgtPrim2Mindm();
		dgtPrim2Mindd = fdgtPrim2Mindd();
		dgtPrim2Mindn = fdgtPrim2Mindn();
		dgtPrim2Maxdm = fdgtPrim2Maxdm();
		dgtPrim2Maxdd = fdgtPrim2Maxdd();
		dgtPrim2Maxdn = fdgtPrim2Maxdn();

		// obliczanie poszczegolnych lambda:
		lbarchiAB = -gt*p.t + (3./4.)*gs*p.J*v.barchiAB;
		lsimchiAB = (1./4.)*gs*p.J*v.simchiAB;
		lchiS = -(1./2.)*p.tPrim*(gtPrimMax + gtPrimMin);
		lchiT = 0.;
		lchiSPrim = -(1./2.)*p.tPrim2*(gtPrim2Max + gtPrim2Min);
		lchiTPrim = 0.;
		lDeltaS = (3./4.)*gs*p.J*v.DeltaS;
		lDeltaT = -(1./4.)*gs*p.J*v.DeltaT;
		lsimDeltaB = -(1./4.)*gs*p.J*v.simDeltaB;
		ln = -16.*p.t*v.barchiAB*dgtdn \
			- 8.*p.tPrim*v.chiS*dgtPrimMindn - 8.*p.tPrim*v.chiS*dgtPrimMaxdn \
			- 8.*p.tPrim2*v.chiSPrim*dgtPrim2Mindn - 8.*p.tPrim2*v.chiSPrim*dgtPrim2Maxdn \
			- 2.*p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB + v.simchiAB*v.simchiAB -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT - 0.5*v.simDeltaB*v.simDeltaB)*dgsdn;
		lmAF = 2.*gs*p.J*v.mAF - 16.*p.t*v.barchiAB*dgtdm \
			- 8.*p.tPrim*v.chiS*dgtPrimMindm - 8.*p.tPrim*v.chiS*dgtPrimMaxdm \
			- 8.*p.tPrim2*v.chiSPrim*dgtPrim2Mindm - 8.*p.tPrim2*v.chiSPrim*dgtPrim2Maxdm \
			- 2.*p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB + v.simchiAB*v.simchiAB -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT - 0.5*v.simDeltaB*v.simDeltaB)*dgsdm;

		// zeby pozniej mozna bylo wypisac
		glob.lbarchiAB = lbarchiAB;
		glob.lsimchiAB = lsimchiAB;
		glob.lchiS = lchiS;
		glob.lchiT = lchiT;
		glob.lchiSPrim = lchiSPrim;
		glob.lchiTPrim = lchiTPrim;
		glob.lDeltaS = lDeltaS;
		glob.lDeltaT = lDeltaT;
		glob.lsimDeltaB = lsimDeltaB;
		glob.ln = ln;
		glob.lmAF = lmAF;
	}

	//REFRESH XY
	void refresh(double kx, double ky) {
		gammak = 2.*(cos(kx) + cos(ky));
		gammaPrimk = 4.*cos(kx)*cos(ky);
		gamma2Primk = 2.*(cos(2.*kx) + cos(2.*ky));
		etak = 2.*(cos(kx) - cos(ky));

		refresh();

		E[0]          		= E1(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlbarchiAB[0]  	= dE1dl_barchiAB(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlsimchiAB[0]  	= dE1dl_simchiAB(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlchiS[0] 		= dE1dl_chiS(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlchiT[0] 		= dE1dl_chiT(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlchiSPrim[0] 	= dE1dl_chiSPrim(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlchiTPrim[0] 	= dE1dl_chiTPrim(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlDeltaS[0] 		= dE1dl_DeltaS(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlDeltaT[0] 		= dE1dl_DeltaT(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlsimDeltaB[0] 	= dE1dl_simDeltaB(gammak, gammaPrimk, gamma2Primk, etak);
		dEdln[0]      		= dE1dl_n(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlmAF[0]    		= dE1dl_mAF(gammak, gammaPrimk, gamma2Primk, etak);

		E[1]          		= -E[0];
		dEdlbarchiAB[1]  	= -dEdlbarchiAB[0];
		dEdlsimchiAB[1]  	= -dEdlsimchiAB[0];
		dEdlchiS[1] 		= -dEdlchiS[0];
		dEdlchiT[1] 		= -dEdlchiT[0];
		dEdlchiSPrim[1] 	= -dEdlchiSPrim[0];
		dEdlchiTPrim[1] 	= -dEdlchiTPrim[0];
		dEdlDeltaS[1] 		= -dEdlDeltaS[0];
		dEdlDeltaT[1] 		= -dEdlDeltaT[0];
		dEdlsimDeltaB[1] 	= -dEdlsimDeltaB[0];
		dEdln[1]      		= -dEdln[0];
		dEdlmAF[1]    		= -dEdlmAF[0];

		E[2]          		= E3(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlbarchiAB[2]  	= dE3dl_barchiAB(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlsimchiAB[2]  	= dE3dl_simchiAB(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlchiS[2] 		= dE3dl_chiS(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlchiT[2] 		= dE3dl_chiT(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlchiSPrim[2] 	= dE3dl_chiSPrim(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlchiTPrim[2] 	= dE3dl_chiTPrim(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlDeltaS[2] 		= dE3dl_DeltaS(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlDeltaT[2] 		= dE3dl_DeltaT(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlsimDeltaB[2] 	= dE3dl_simDeltaB(gammak, gammaPrimk, gamma2Primk, etak);
		dEdln[2]      		= dE3dl_n(gammak, gammaPrimk, gamma2Primk, etak);
		dEdlmAF[2]    		= dE3dl_mAF(gammak, gammaPrimk, gamma2Primk, etak);

		E[3]          		= -E[2];
		dEdlbarchiAB[3]  	= -dEdlbarchiAB[2];
		dEdlsimchiAB[3]  	= -dEdlsimchiAB[2];
		dEdlchiS[3] 		= -dEdlchiS[2];
		dEdlchiT[3] 		= -dEdlchiT[2];
		dEdlchiSPrim[3] 	= -dEdlchiSPrim[2];
		dEdlchiTPrim[3] 	= -dEdlchiTPrim[2];
		dEdlDeltaS[3] 		= -dEdlDeltaS[2];
		dEdlDeltaT[3] 		= -dEdlDeltaT[2];
		dEdlsimDeltaB[3] 	= -dEdlsimDeltaB[2];
		dEdln[3]      		= -dEdln[2];
		dEdlmAF[3]    		= -dEdlmAF[2];
	}
};


/////////////////
///////////////// FullSimplify
/////////////////

double Func::fgBDelta() {
    //n = 1.-p.sdelta;
    double r = n/2. + v.mAF/2.;
    double w = n/2. - v.mAF/2.;
    double d=v.d*v.d;
	double part1 = (n-2.*d)/(2.*(n-2.*r*w));
	double part2 = sqrt((1.-w)*(1.-n+d)/(1.-r)) + sqrt(w*d/r);
	double part3 = sqrt((1.-r)*(1.-n+d)/(1.-w)) + sqrt(r*d/w);
	return part1*(part2*part2 + part3*part3);
}

double Func::fgm() {
   double r = n/2. + v.mAF/2.;
   double w = n/2. - v.mAF/2.;
   double d = v.d*v.d;
   return (n-2.*d)/(n-2.*w*r);
}

double Func::fgt() {
	return (2*(2*pow(d,2) - n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))))/
			   (-pow(m,2) + (-2 + n)*n);
}

double Func::fgs() {
	return (4*pow(-2*pow(d,2) + n,2))/pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fgtPrimMin() {
	return (2*(2*pow(d,2) - n)*pow(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)),2))/(-pow(m,2) + (-2 + n)*n);
}

double Func::fgtPrimMax() {
	return (2*(2*pow(d,2) - n)*pow(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)),2))/(-pow(m,2) + (-2 + n)*n);
}

double Func::fgtPrim2Min() {
	return fgtPrimMin();
}

double Func::fgtPrim2Max() {
	return fgtPrimMax();
}

////////////////////////////////

double Func::fdgtdm() {
	return (4*(-2*pow(d,2) + n)*(-(m*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
	          (sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)))) +
	       (pow(m,4)*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n))*sqrt((pow(d,2)*(-m + n))/(m + n)) - sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n))*sqrt((pow(d,2)*(m + n))/(-m + n))))/
	        ((m - n)*(2 + m - n)*(-2 + m + n)*(m + n)) - (pow(-2 + n,2)*pow(n,2)*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n))*sqrt((pow(d,2)*(-m + n))/(m + n)) -
	            sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n))*sqrt((pow(d,2)*(m + n))/(-m + n))))/((m - n)*(2 + m - n)*(-2 + m + n)*(m + n))))/pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgtdd() {
	return (2*(-4*d*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))) +
		       (-2*pow(d,2) + n)*((d*sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)))/(1 + pow(d,2) - n) + sqrt((pow(d,2)*(-m + n))/(m + n))/d)*
		        (sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))) +
		       (-2*pow(d,2) + n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
		        ((d*sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)))/(1 + pow(d,2) - n) + sqrt((pow(d,2)*(m + n))/(-m + n))/d)))/(pow(m,2) - (-2 + n)*n);
}

double Func::fdgtdn() {
	return (4*(((pow(m,2) - (-2 + n)*n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) +
            sqrt((pow(d,2)*(-m + n))/(m + n)))*
          (sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt(-((pow(d,2)*(m + n))/(m - n)))))/2.
         - (1 - n)*(-2*pow(d,2) + n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) +
          sqrt((pow(d,2)*(-m + n))/(m + n)))*
        (sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))) +
       ((-2*pow(d,2) + n)*(pow(m,2) - (-2 + n)*n)*
          ((pow(m,2) + 2*m*(1 + pow(d,2) - n) - pow(-2 + n,2))/
             (sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n))*pow(-2 + m + n,2)) +
            (2*pow(d,2)*m)/(sqrt((pow(d,2)*(-m + n))/(m + n))*pow(m + n,2)))*
          (sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))))/4.\
        + ((-2*pow(d,2) + n)*(pow(m,2) - (-2 + n)*n)*
          (sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
          ((pow(m,2) - 2*m*(1 + pow(d,2) - n) - pow(-2 + n,2))/
             (pow(2 + m - n,2)*sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n))) +
            (2*m*sqrt(-((pow(d,2)*(m + n))/(m - n))))/(pow(m,2) - pow(n,2))))/4.))/
   pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgsdm() {
	return (-16*m*pow(-2*pow(d,2) + n,2))/pow(pow(m,2) - (-2 + n)*n,3);
}

double Func::fdgsdd() {
	return (64*pow(d,3) - 32*d*n)/pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgsdn() {
	return (8*(-2*pow(d,2) + n)*(pow(m,2) - 4*pow(d,2)*(-1 + n) + pow(n,2)))/pow(pow(m,2) - (-2 + n)*n,3);
}

double Func::fdgtPrimMindm() {
	return (4*(-2*pow(d,2) + n)*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)))*
		     (-(m*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)))) +
		       ((pow(m,2) - (-2 + n)*n)*((2*(-2 + n)*sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)))/((-2 - m + n)*(-2 + m + n)) + (2*n*sqrt((pow(d,2)*(m + n))/(-m + n)))/(-pow(m,2) + pow(n,2))))/2.))/
		   pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgtPrimMindd() {
	return (4*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)))*
		     (-2*d*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))) +
		       (-2*pow(d,2) + n)*((d*sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)))/(1 + pow(d,2) - n) + sqrt((pow(d,2)*(m + n))/(-m + n))/d)))/(pow(m,2) - (-2 + n)*n);
}

double Func::fdgtPrimMindn() {
	return (4*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)))*
		     (((pow(m,2) - (-2 + n)*n)*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt(-((pow(d,2)*(m + n))/(m - n)))))/2. +
		       (-1 + n)*(-2*pow(d,2) + n)*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))) +
		       ((-2*pow(d,2) + n)*(pow(m,2) - (-2 + n)*n)*((pow(m,2) - 2*m*(1 + pow(d,2) - n) - pow(-2 + n,2))/(pow(2 + m - n,2)*sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n))) +
		            (2*m*sqrt(-((pow(d,2)*(m + n))/(m - n))))/(pow(m,2) - pow(n,2))))/2.))/pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgtPrimMaxdm() {
	return (4*(-2*pow(d,2) + n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
		     (-(m*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))) +
		       ((pow(m,2) - (-2 + n)*n)*((2*(-2 + n)*sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)))/((2 + m - n)*(-2 + m + n)) - (2*pow(d,2)*n)/(sqrt((pow(d,2)*(-m + n))/(m + n))*pow(m + n,2))))/2.))/
		   pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgtPrimMaxdd() {
	return (4*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
		     (-2*d*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n))) +
		       (-2*pow(d,2) + n)*((d*sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)))/(1 + pow(d,2) - n) + sqrt((pow(d,2)*(-m + n))/(m + n))/d)))/(pow(m,2) - (-2 + n)*n);
}

double Func::fdgtPrimMaxdn() {
	return (4*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
		     ((-1 + n)*(-2*pow(d,2) + n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n))) +
		       ((pow(m,2) - (-2 + n)*n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n))))/2. +
		       ((-2*pow(d,2) + n)*(pow(m,2) - (-2 + n)*n)*((pow(m,2) + 2*m*(1 + pow(d,2) - n) - pow(-2 + n,2))/(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n))*pow(-2 + m + n,2)) +
		            (2*pow(d,2)*m)/(sqrt((pow(d,2)*(-m + n))/(m + n))*pow(m + n,2))))/2.))/pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgtPrim2Mindm() {
	return fdgtPrimMindm();
}

double Func::fdgtPrim2Mindd() {
	return fdgtPrimMindd();
}

double Func::fdgtPrim2Mindn() {
	return fdgtPrimMindn();
}

double Func::fdgtPrim2Maxdm() {
	return fdgtPrimMaxdm();
}

double Func::fdgtPrim2Maxdd() {
	return fdgtPrimMaxdd();
}

double Func::fdgtPrim2Maxdn() {
	return fdgtPrimMaxdn();
}

/////////////////////////////////////////////////////////

double Func::E1(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) + 4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) +
		      2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) + pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
		      4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		           (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		         4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*(pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))/2.;
}

double Func::E3(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) + 4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) +
		      pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
		      2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
		         sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) + (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		           4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*(pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/2.;
}

double Func::dE1dl_barchiAB(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return (2*pow(gammak,2)*lbarchiAB*(pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2) -
		       sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		           (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		         4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		          (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
		   (sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		         (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		        (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))*
		     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + 2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
		       pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
		       4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))));
}

double Func::dE3dl_barchiAB(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return (-2*pow(gammak,2)*lbarchiAB*(pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2) +
		       sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		           (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		         4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		          (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
		   (sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		         (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		        (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))*
		     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(ln + 2*mu,2) +
		       4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
		       2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
		          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))));
}

double Func::dE1dl_simchiAB(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return (2*pow(gammak,2)*lsimchiAB*(pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2) -
		       sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		           (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		         4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		          (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
		   (sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		         (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		        (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))*
		     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + 2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
		       pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
		       4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))));
}

double Func::dE3dl_simchiAB(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return (-2*pow(gammak,2)*lsimchiAB*(pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2) +
		       sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		           (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		         4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		          (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
		   (sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		         (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		        (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))*
		     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(ln + 2*mu,2) +
		       4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
		       2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
		          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))));
}

double Func::dE1dl_chiS(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((gammaPrimk*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu +
	         (-4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu) -
	            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	               (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + 2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	       pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
	       4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))));
}

double Func::dE3dl_chiS(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((gammaPrimk*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu +
	         (4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu) +
	            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	               (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(ln + 2*mu,2) +
	       4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
	       2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))));
}

double Func::dE1dl_chiT(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((gammaPrimk*(2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF -
	         ((2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)*
	            (pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + 2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	       pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
	       4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))));
}

double Func::dE3dl_chiT(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((gammaPrimk*(2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF +
	         ((2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)*
	            (pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(ln + 2*mu,2) +
	       4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
	       2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))));
}

double Func::dE1dl_chiSPrim(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((gamma2Primk*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu +
	         (-4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu) -
	            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	               (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + 2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	       pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
	       4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))));
}

double Func::dE3dl_chiSPrim(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((gamma2Primk*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu +
	         (4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu) +
	            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	               (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(ln + 2*mu,2) +
	       4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
	       2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))));
}

double Func::dE1dl_chiTPrim(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((gamma2Primk*(2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF -
	         ((2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)*
	            (pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + 2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	       pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
	       4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))));
}

double Func::dE3dl_chiTPrim(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((gamma2Primk*(2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF +
	         ((2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)*
	            (pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(ln + 2*mu,2) +
	       4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
	       2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))));
}

double Func::dE1dl_DeltaS(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return (2*pow(etak,2)*(pow(etak,2)*lDeltaT*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		       lDeltaT*(2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu) -
		       lDeltaS*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
		   (sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		         (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		        (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))*
		     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + 2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
		       pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
		       4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))));
}

double Func::dE3dl_DeltaS(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return (2*pow(etak,2)*(pow(etak,2)*lDeltaT*(-4*lDeltaS*lDeltaT + pow(lsimDeltaB,2)) -
		       lDeltaT*(2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu) -
		       lDeltaS*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
		   (sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		         (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		        (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))*
		     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(ln + 2*mu,2) +
		       4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
		       2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
		          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))));
}

double Func::dE1dl_DeltaT(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return (2*pow(etak,2)*(4*pow(gammaPrimk,2)*lchiS*lchiT*lDeltaS + 4*pow(gamma2Primk,2)*lchiSPrim*lchiTPrim*lDeltaS + 4*pow(gammak,2)*pow(lbarchiAB,2)*lDeltaT +
		       4*pow(etak,2)*pow(lDeltaS,2)*lDeltaT + lDeltaS*lmAF*ln + 4*pow(gammak,2)*lDeltaT*pow(lsimchiAB,2) - pow(etak,2)*lDeltaS*pow(lsimDeltaB,2) +
		       2*lDeltaS*lmAF*mu + 2*gammaPrimk*lDeltaS*(2*gamma2Primk*(lchiSPrim*lchiT + lchiS*lchiTPrim) + lchiS*lmAF + lchiT*ln + 2*lchiT*mu) +
		       2*gamma2Primk*lDeltaS*(lchiSPrim*lmAF + lchiTPrim*(ln + 2*mu)) -
		       lDeltaT*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
		   (sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		         (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		        (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))*
		     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + 2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
		       pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
		       4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))));
}

double Func::dE3dl_DeltaT(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return (-2*pow(etak,2)*(4*pow(gammaPrimk,2)*lchiS*lchiT*lDeltaS + 4*pow(gamma2Primk,2)*lchiSPrim*lchiTPrim*lDeltaS + 4*pow(gammak,2)*pow(lbarchiAB,2)*lDeltaT +
		       4*pow(etak,2)*pow(lDeltaS,2)*lDeltaT + lDeltaS*lmAF*ln + 4*pow(gammak,2)*lDeltaT*pow(lsimchiAB,2) - pow(etak,2)*lDeltaS*pow(lsimDeltaB,2) +
		       2*lDeltaS*lmAF*mu + 2*gammaPrimk*lDeltaS*(2*gamma2Primk*(lchiSPrim*lchiT + lchiS*lchiTPrim) + lchiS*lmAF + lchiT*ln + 2*lchiT*mu) +
		       2*gamma2Primk*lDeltaS*(lchiSPrim*lmAF + lchiTPrim*(ln + 2*mu)) +
		       lDeltaT*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
		   (sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		         (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		        (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))*
		     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
		       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(ln + 2*mu,2) +
		       4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
		       2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
		          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
		              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
		            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
		             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))));
}

double Func::dE1dl_simDeltaB(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((pow(etak,2)*lsimDeltaB*(-2*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	         (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu) +
	         sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	             (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	           4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	            (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     (sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	           (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	         4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	          (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))*
	       sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	         4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + 2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	         pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
	         4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))));
}

double Func::dE3dl_simDeltaB(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((pow(etak,2)*lsimDeltaB*(2*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(etak,2)*(-4*lDeltaS*lDeltaT + pow(lsimDeltaB,2)) -
	         (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu) +
	         sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	             (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	           4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	            (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     (sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	           (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	         4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	          (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))*
	       sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	         4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(ln + 2*mu,2) +
	         4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
	         2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	            sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	                (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	              4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	               (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))));
}

double Func::dE1dl_n(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((gammaPrimk*lchiS + gamma2Primk*lchiSPrim + ln/2. + mu + (-4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	           (2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu) -
	          (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	             (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	        (2.*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + 2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	       pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
	       4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))));
}

double Func::dE3dl_n(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((gammaPrimk*lchiS + gamma2Primk*lchiSPrim + ln/2. + mu + (4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	           (2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu) +
	          (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	             (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	        (2.*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(ln + 2*mu,2) +
	       4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
	       2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))));
}

double Func::dE1dl_mAF(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -((gammaPrimk*lchiT + gamma2Primk*lchiTPrim + lmAF/2. - ((2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)*
	          (pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	        (2.*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))))/
	     sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + 2*pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	       pow(ln + 2*mu,2) + 4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) +
	       4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) - 2*sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	            (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	          4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	           (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2)))));
}

double Func::dE3dl_mAF(double gammak, double gammaPrimk, double gamma2Primk, double etak) {
	return -(2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF + ((2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)*
	         (pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	           (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu)))/
	       sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	           (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	         4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	          (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))/
	   (2.*sqrt(4*pow(gammaPrimk,2)*(pow(lchiS,2) + pow(lchiT,2)) + 4*pow(gamma2Primk,2)*(pow(lchiSPrim,2) + pow(lchiTPrim,2)) + pow(lmAF,2) +
	       4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2)) + pow(ln + 2*mu,2) +
	       4*gammaPrimk*(2*gamma2Primk*(lchiS*lchiSPrim + lchiT*lchiTPrim) + lchiT*lmAF + lchiS*(ln + 2*mu)) + 4*gamma2Primk*(lchiTPrim*lmAF + lchiSPrim*(ln + 2*mu)) +
	       2*(pow(etak,2)*(2*(pow(lDeltaS,2) + pow(lDeltaT,2)) + pow(lsimDeltaB,2)) +
	          sqrt(pow(pow(etak,2)*(4*lDeltaS*lDeltaT - pow(lsimDeltaB,2)) +
	              (2*gammaPrimk*lchiT + 2*gamma2Primk*lchiTPrim + lmAF)*(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu),2) +
	            4*pow(gammak,2)*(pow(lbarchiAB,2) + pow(lsimchiAB,2))*
	             (pow(etak,2)*(4*pow(lDeltaT,2) + pow(lsimDeltaB,2)) + pow(2*gammaPrimk*lchiS + 2*gamma2Primk*lchiSPrim + ln + 2*mu,2))))));
}


/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	void onePoint(variables &y, int i, int j);

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->set(p);
	}

	void run(variables &w);
};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	f->refresh(x, y);
	if(DEEP_DEBUG) { f->printLambda(); }

	double Etan[4] = { (1-tanh(p.beta*f->E[0]/2.))/2., (1-tanh(p.beta*f->E[1]/2.))/2., (1-tanh(p.beta*f->E[2]/2.))/2., (1-tanh(p.beta*f->E[3]/2.))/2. };

	for(int i=0; i<4; i++) {
		w.barchiAB += Etan[i]*f->dEdlbarchiAB[i];
		w.simchiAB += Etan[i]*f->dEdlsimchiAB[i];
		w.chiS += Etan[i]*f->dEdlchiS[i];
		w.chiT += Etan[i]*f->dEdlchiT[i];
		w.chiSPrim += Etan[i]*f->dEdlchiSPrim[i];
		w.chiTPrim += Etan[i]*f->dEdlchiTPrim[i];
		w.DeltaS += Etan[i]*f->dEdlDeltaS[i];
		w.DeltaT += Etan[i]*f->dEdlDeltaT[i];
		w.simDeltaB += Etan[i]*f->dEdlsimDeltaB[i];
		w.mu += Etan[i]*f->dEdln[i];
		w.mAF += Etan[i]*f->dEdlmAF[i];
		// w.d w tym miejscu jeszcze nie

        //if(f->E[i]>0) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
        //if(f->E[i]<0) { w.Energy += f->E[i]; }

        // it would be safer have following code - warning: 30 is chosen arbitrarily
		if(p.beta*f->E[i] > 30.) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
		else if(p.beta*f->E[i] < -30) { w.Energy += f->E[i]; }
		else { w.Energy += -log(1 + exp(-p.beta*f->E[i]))/p.beta; }
	}
}

void Equa::run(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation
	unsigned int max = p.L/2; // attention: L should be even!

	// SERIAL CASE
	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
	// * edges (2 times)
	for(unsigned int i=1; i<max; i++) {
		onePoint(temp, i, 0);
		onePoint(temp, 0, i);
		onePoint(temp, i, max-i);
	}
   addTo(temp, w, 2);
   clear(temp);
   // * interior area (4 times) - SERIAL CASE
   for(unsigned int i=1; i<max; i++) {
		for(unsigned int j=1; j<i; j++) {
			onePoint(temp, max-i, j);
     	}
	}
	addTo(temp, w, 4);
	clear(temp);
	// * missing corners
	onePoint(w, max, 0);
	//onePoint(w, 0, max);
	onePoint(w, 0, 0);

	// what was out of sumation
	double L2 = (double)(p.L*p.L);
	w.barchiAB /= L2;
	w.simchiAB /= L2;
	w.chiS /= L2;
	w.chiT /= L2;
	w.chiSPrim /= L2;
	w.chiTPrim /= L2;
	w.DeltaS /= L2;
	w.DeltaT /= L2;
	w.simDeltaB /= L2;
	w.mu /= L2;
	w.mAF /= L2;
    w.Energy /= L2;

	w.barchiAB += 8.*v.barchiAB;
	w.simchiAB += -8.*v.simchiAB;
	w.chiS += 8.*v.chiS;
	w.chiT += 8.*v.chiT;
	w.chiSPrim += 8.*v.chiSPrim;
	w.chiTPrim += 8.*v.chiTPrim;
	w.DeltaS += 8.*v.DeltaS;
	w.DeltaT += 8.*v.DeltaT;
	w.simDeltaB += -4.*v.simDeltaB;
	w.mu += -0.5*(1-p.n);
	w.mAF += 0.5*v.mAF;

	// calculate dd
	w.d = 2.*v.d*p.U \
			+ 8.*p.t*v.barchiAB*f->dgtdd \
			+ 4.*p.tPrim*v.chiS*f->dgtPrimMaxdd + 4.*p.tPrim*v.chiS*f->dgtPrimMindd \
			+ 4.*p.tPrim2*v.chiSPrim*f->dgtPrim2Maxdd + 4.*p.tPrim2*v.chiSPrim*f->dgtPrim2Mindd \
			+ p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB + v.simchiAB*v.simchiAB -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT - 0.5*v.simDeltaB*v.simDeltaB)*f->dgsdd;

    v.Energy = w.Energy \
    		+ 8.*f->gt*p.t*v.barchiAB \
    		+ 4.*f->gtPrimMax*p.tPrim*v.chiS + 4.*f->gtPrimMin*p.tPrim*v.chiS \
    		+ 4.*f->gtPrim2Max*p.tPrim2*v.chiSPrim + 4.*f->gtPrim2Min*p.tPrim2*v.chiSPrim \
    		+ f->gs*p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB + v.simchiAB*v.simchiAB -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT - 0.5*v.simDeltaB*v.simDeltaB) \
    		+ 8.*(f->lbarchiAB*v.barchiAB - f->lsimchiAB*v.simchiAB + f->lchiS*v.chiS + f->lchiT*v.chiT + f->lchiSPrim*v.chiSPrim + f->lchiTPrim*v.chiTPrim +\
    				f->lDeltaS*v.DeltaS + f->lDeltaT*v.DeltaT - 0.5*f->lsimDeltaB*v.simDeltaB) \
    		+ p.U*v.d*v.d + 0.5*f->ln*(p.n-1) + 0.5*f->lmAF*v.mAF -v.mu;

	//if(w.chiAB != w.chiAB) { w.chiAB = 1000; }
	//if(w.chiS != w.chiS) { w.chiS = 1000; }
	//if(w.chiT != w.chiT) { w.chiT = 1000; }
	//if(w.barDelta != w.barDelta) { w.barDelta = 1000; }
	//if(w.simDelta != w.simDelta) { w.simDelta = 1000; }
	//if(w.mu != w.mu) { w.mu = 1000; }
	//if(w.mAF != w.mAF) { w.mAF = 1000; }

    // results
    case_2D3t_SGA::glob.Energy = v.Energy + v.mu*p.n;
    case_2D3t_SGA::glob.gt = f->fgt();
    case_2D3t_SGA::glob.gs = f->fgs();
    case_2D3t_SGA::glob.gtPrimMax = f->fgtPrimMax();
    case_2D3t_SGA::glob.gtPrimMin = f->fgtPrimMin();
    case_2D3t_SGA::glob.gtPrim2Max = f->fgtPrim2Max();
    case_2D3t_SGA::glob.gtPrim2Min = f->fgtPrim2Min();
	case_2D3t_SGA::glob.gs = f->fgs();
	case_2D3t_SGA::glob.gBDelta = f->fgBDelta();
	case_2D3t_SGA::glob.gm = f->fgm();
}

// ***********************************************************************************  //
// *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY ***
// ***********************************************************************************  //


int function(const gsl_vector * x, void *params, gsl_vector * f) {

	struct parameters p;
	p.n = ((Par *) params)->n();
	p.J = ((Par *) params)->J();
	p.t = ((Par *) params)->t();
	p.tPrim = ((Par *) params)->tPrim();
	p.tPrim2 = ((Par *) params)->tPrim2();
	p.U = ((Par *) params)->U();
	p.beta = ((Par *) params)->beta();
	p.L = ((Par *) params)->L();

	variables v;
	v.barchiAB = gsl_vector_get (x, 0);
	v.simchiAB = gsl_vector_get (x, 1);
	v.chiS = gsl_vector_get (x, 2);
	v.chiT = gsl_vector_get (x, 3);
	v.chiSPrim = gsl_vector_get (x, 4);
	v.chiTPrim = gsl_vector_get (x, 5);
	v.DeltaS = gsl_vector_get (x, 6);
	v.DeltaT = gsl_vector_get (x, 7);
	v.simDeltaB = gsl_vector_get (x, 8);
	v.mu = gsl_vector_get (x, 9);
	v.mAF = gsl_vector_get (x, 10);
	v.d = gsl_vector_get (x, 11);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7), gsl_vector_get (x, 8),
				gsl_vector_get (x, 9), gsl_vector_get (x, 10), gsl_vector_get (x, 11));
	}

    // It's start of the calculations

    Equa E(p, v);
    variables w;
    E.run(w);

    // It's end of the calculations

	gsl_vector_set (f, 0, w.barchiAB);
	gsl_vector_set (f, 1, w.simchiAB);
	gsl_vector_set (f, 2, w.chiS);
	gsl_vector_set (f, 3, w.chiT);
	gsl_vector_set (f, 4, w.chiSPrim);
	gsl_vector_set (f, 5, w.chiTPrim);
	gsl_vector_set (f, 6, w.DeltaS);
	gsl_vector_set (f, 7, w.DeltaT);
	gsl_vector_set (f, 8, w.simDeltaB);
	gsl_vector_set (f, 9, w.mu);
	gsl_vector_set (f, 10, w.mAF);
	gsl_vector_set (f, 11, w.d);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7), gsl_vector_get (x, 8),
				gsl_vector_get (x, 9), gsl_vector_get (x, 10), gsl_vector_get (x, 11));
		printf("###  %24.15e# %25.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# | E = %24.15e\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3), gsl_vector_get (f, 4), gsl_vector_get (f, 5),
			gsl_vector_get (f, 6), gsl_vector_get (f, 7), gsl_vector_get (f, 8),
			gsl_vector_get (f, 9), gsl_vector_get (f, 10), gsl_vector_get (f, 11), case_2D3t_SGA::glob.Energy);
		printf("\n");
	}

    return GSL_SUCCESS;
}

// ************* CLASSES ************* //

// Par

void Par::print(FILE *file) const {
	for(int i=0; i<NtabD; i++) { fprintf(file, " %23.15e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
}

void Par::print_verboten(FILE *file) const {
	fprintf(file, "PARAMETERS: n=%f, J=%f, t=%f, t'=%f, t''=%f, U=%f, beta=%f, L=%d\n", n(), J(), t(), tPrim(), tPrim2(), U(), beta(), L());
}

void Par::print_label(FILE *file) const {
	fprintf(file, "n   J   t   t'   t''   U   beta   L\n");
}

// Var

void Var::print(FILE *file) const{
	for(int i=0; i<NtabD; i++) { fprintf(file, " %23.15e", tabD[i]); }
}

void Var::print_verboten(FILE *file) const{
	fprintf(file, "VARIABLES:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %23.15e", tabD[i]); }
	fprintf(file, "\n");
}

void Var::print_label(FILE *file) const {
	fprintf(file, "    barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d");
}

// Res

void Res::set(gsl_multiroot_fsolver *s) {
	for(int i=0; i<n; i++) {
		(*V)[i] = gsl_vector_get (s->x, i);
	}
	this->glob = case_2D3t_SGA::glob;
}

void Res::print(FILE *file) const { V->print(file); fprintf(file, "%23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f\n",
		glob.Energy, glob.gt, glob.gtPrimMax, glob.gtPrimMin, glob.gtPrim2Max, glob.gtPrim2Min, glob.gs, glob.gBDelta, glob.gm, glob.lbarchiAB, glob.lsimchiAB, glob.lchiS, glob.lchiT, glob.lchiSPrim, glob.lchiTPrim, glob.lDeltaS, glob.lDeltaT, glob.lsimDeltaB, glob.ln, glob.lmAF); }
void Res::print_verboten(FILE *file) const { fprintf(file, "Result: "); V->print_verboten(file); fprintf(file, "%23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f %23.15f\n", glob.Energy, glob.gt, glob.gtPrimMax, glob.gtPrimMin, glob.gtPrim2Max, glob.gtPrim2Min, glob.gs, glob.gBDelta, glob.gm, glob.lbarchiAB, glob.lsimchiAB, glob.lchiS, glob.lchiT, glob.lchiSPrim, glob.lchiTPrim, glob.lDeltaS, glob.lDeltaT, glob.lsimDeltaB, glob.ln, glob.lmAF); }
void Res::print_label(FILE *file) const { V->print_label(file); fprintf(file, "                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF\n");}

void example() {

	printf("\nEXAMPLE 2D2t_fast\n");

	Par par(N_PAR);
	par.set_n(0.999);
	par.set_J(1./3.);
	par.set_t(-1.);
	par.set_tPrim(0.);
	par.set_tPrim2(-0.);
	par.set_U(12);
	par.set_beta(1500);
	par.set_L(256);
	printf("#"); par.print_verboten(); fflush(stdout);

    double ZERO = 1e-12;

    Var initV_Z2_U12_n999_B1500(N_VAR);
    initV_Z2_U12_n999_B1500.set_barchiAB(1.693217473268961e-01);
    initV_Z2_U12_n999_B1500.set_simchiAB(ZERO);
    initV_Z2_U12_n999_B1500.set_chiS(1.442308124273939e-04);
    initV_Z2_U12_n999_B1500.set_chiT(-3.865254347226006e-04);
    initV_Z2_U12_n999_B1500.set_chiSPrim(8.791624100731755e-05);
    initV_Z2_U12_n999_B1500.set_chiTPrim(-1.434708493837335e-02);
    initV_Z2_U12_n999_B1500.set_DeltaS(1.669293228458277e-01);
    initV_Z2_U12_n999_B1500.set_DeltaT(3.915289953743466e-05);
    initV_Z2_U12_n999_B1500.set_simDeltaB(ZERO);
    initV_Z2_U12_n999_B1500.set_mu(3.371567247648004e+00);
    initV_Z2_U12_n999_B1500.set_mAF(1.317210700617173e-01);
    initV_Z2_U12_n999_B1500.set_d(1.496803286385485e-02);

    printf("Starting point: "); initV_Z2_U12_n999_B1500.print(); printf("\n\n");
    Res res(N_VAR); res.print_label();
    RootGSL root(par, function, res);
    if(root.Solve(initV_Z2_U12_n999_B1500) == GSL_SUCCESS) { root.print(); }

    printf("Koniec\n");
}

void example_1D_beta_run(Par &par, Var &initV) {
	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	initV.shake();
	if(root.Solve(initV) != GSL_SUCCESS) { printf("#D# First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n           "); res.print_label();
	fflush(stdout);

	int MNOZNIK = 1;
	/*	LOOPS = 6*MNOZNIK;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_beta(1500 - (double)i*200/MNOZNIK);
		printf("%f ", par.beta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n");

	MNOZNIK = 1;
	LOOPS = 10*MNOZNIK;
	for(int i=0; i<=LOOPS; i++) {
		par.set_beta(300 - (double)i*20/MNOZNIK);
		printf("%f ", par.beta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n");

	MNOZNIK = 1;
	LOOPS = 12*MNOZNIK;
	for(int i=0; i<=LOOPS; i++) {
		par.set_beta(100 - (double)i*5/MNOZNIK);
		printf("%f ", par.beta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n");
*/
	MNOZNIK = 1;
	LOOPS = 50*MNOZNIK;
	for(int i=0; i<=LOOPS; i++) {
		par.set_beta(50. - (double)i*1/MNOZNIK);
		printf("%f ", par.beta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
}

void example_1D_sdelta_run(Par &par, Var &initV) {
	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	initV.shake();
	if(root.Solve(initV) != GSL_SUCCESS) { printf("#D# First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n           "); res.print_label();
	fflush(stdout);

	int MNOZNIK = 1;
	LOOPS = 10*MNOZNIK;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.999 + (double)i*0.0001/MNOZNIK);
		printf("%f ", par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 90*MNOZNIK;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.999 - (double)i*0.0001/MNOZNIK);
		printf("%f ", par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only resolution. Still the same case\n");

	LOOPS = 590*MNOZNIK;
	for(int i=0; i<=LOOPS; i++) {
		par.set_n(0.99 - (double)i*0.001/MNOZNIK);
		printf("%f ", par.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
}

void example_1D() {
    Par par(N_PAR);
    par.set_n(0.999);
    par.set_J(1./3.);
    par.set_t(-1.);
    par.set_tPrim(0.25);
    par.set_tPrim2(0.);
    par.set_U(12);
    par.set_beta(10.);
    par.set_L(256);
    par.print_verboten();

    double ZERO = 1e-15;

//    PARAMETERS: n=0.998000, J=0.333333, t=-1.000000, t'=0.000000, t''=0.000000, U=12.000000, beta=500.000000, L=256
//    INTERESUJĄCE
//    sdelta                    barchiAB             simchiAB                 chiS                 chiT                     chiSPrim                chiTPrim               DeltaS                  DeltaT                simDeltaB               mu                     mAF                    d                      E                      gt               gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                      gBDelta               gm                 lbarchiAB                 lsimchiAB             lDeltaS              lchiS              lchiT                  lchiSPrim              lchiTPrim              lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
//    DS                  1.715797881677855e-01 -1.126341666310388e-23 2.672988098584522e-04 -2.437895155902840e-14  2.477049655097007e-04 -4.472539643978754e-13 -1.670753025582039e-01 -7.648477242122892e-14  2.546951936974829e-20 3.460201913896650e+00  3.954157982369582e-12  2.168172760566418e-02   -0.236375951167834     0.010160801009590     0.010160801009622     0.010160801009559     0.010160801009622     0.010160801009559     3.976544834281297     0.010160801009590     1.994127587262484     0.180734481086012    -0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000     0.000000000000000    -0.166095607830953     0.000000000000025    -0.000000000000000    -6.923475012518667     0.000000000003234
//    DS+DT+chiT+m        1.681305657107867e-01  1.731831773264103e-21 6.107227423471418e-04 -4.861621154282696e-02 -3.244360647922483e-04 -1.598368940936969e-02  4.672072595867999e-02  2.175818210965474e-03  5.005692372859998e-14 5.391034757927860e+00  6.509518143510621e-01  8.233828080112927e-02   -0.211788412575583     0.074236893593750     0.080384036761959     0.068559835913315     0.080384036761959     0.068559835913315     1.912411468857998     0.074471936337637     1.382899659721557     0.154620599126473     0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000     0.000000000000000     0.022337313039188    -0.000346754975067    -0.000000000000008   -10.813429193616289     0.390410667317855
//    0                   2.026312528393326e-01  4.227990941756281e-25 6.891768026409183e-04  8.628367873787484e-14 -3.783852286645283e-04 -4.073204472344117e-14 -3.552596135747061e-18  6.609526106245610e-17  2.857903213852130e-13 5.563039518274006e+00 -4.228220586641046e-13  7.975516757789704e-02   -0.176758926698306     0.057751449458077     0.057751449458074     0.057751449458080     0.057751449458074     0.057751449458080     3.883123799303418     0.057751449458077     1.970564335235827     0.254462009553847     0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000    -0.000000000000000    -0.000000000000092   -11.127243656127373    -0.000000000000555

    // I'm choosing start point
    Var initV_AA(N_VAR);			// DS
//    sdelta                    barchiAB             simchiAB                 chiS                 chiT                     chiSPrim                chiTPrim               DeltaS                  DeltaT                simDeltaB               mu                     mAF                    d                      E                      gt               gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                      gBDelta               gm                 lbarchiAB                 lsimchiAB             lDeltaS              lchiS              lchiT                  lchiSPrim              lchiTPrim              lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
//    DS                  1.715797881677855e-01 -1.126341666310388e-23 2.672988098584522e-04 -2.437895155902840e-14  2.477049655097007e-04 -4.472539643978754e-13 -1.670753025582039e-01 -7.648477242122892e-14  2.546951936974829e-20 3.460201913896650e+00  3.954157982369582e-12  2.168172760566418e-02   -0.236375951167834     0.010160801009590     0.010160801009622     0.010160801009559     0.010160801009622     0.010160801009559     3.976544834281297     0.010160801009590     1.994127587262484     0.180734481086012    -0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000     0.000000000000000    -0.166095607830953     0.000000000000025    -0.000000000000000    -6.923475012518667     0.000000000003234
    // DS
    initV_AA.set_barchiAB(1.715797881677855e-01);
    initV_AA.set_simchiAB(ZERO);
    initV_AA.set_chiS(2.672988098584522e-04);
    initV_AA.set_chiT(ZERO);
    initV_AA.set_chiSPrim(2.477049655097007e-04);
    initV_AA.set_chiTPrim(ZERO);
    initV_AA.set_DeltaS(1.670753025582039e-01);
    initV_AA.set_DeltaT(ZERO);
    initV_AA.set_simDeltaB(ZERO);
    initV_AA.set_mu(3.460201913896650e+00);
    initV_AA.set_mAF(ZERO);
    initV_AA.set_d(2.168172760566418e-02);

//    sdelta                    barchiAB             simchiAB                 chiS                 chiT                     chiSPrim                chiTPrim               DeltaS                  DeltaT                simDeltaB               mu                     mAF                    d                      E                      gt               gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                      gBDelta               gm                 lbarchiAB                 lsimchiAB             lDeltaS              lchiS              lchiT                  lchiSPrim              lchiTPrim              lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
//    DS+DT+chiT+m        1.681305657107867e-01  1.731831773264103e-21 6.107227423471418e-04 -4.861621154282696e-02 -3.244360647922483e-04 -1.598368940936969e-02  4.672072595867999e-02  2.175818210965474e-03  5.005692372859998e-14 5.391034757927860e+00  6.509518143510621e-01  8.233828080112927e-02   -0.211788412575583     0.074236893593750     0.080384036761959     0.068559835913315     0.080384036761959     0.068559835913315     1.912411468857998     0.074471936337637     1.382899659721557     0.154620599126473     0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000     0.000000000000000     0.022337313039188    -0.000346754975067    -0.000000000000008   -10.813429193616289     0.390410667317855
    // DS+DT+chiT+m
    Var initV_BB(N_VAR);
    initV_BB.set_barchiAB(1.681305657107867e-01);
    initV_BB.set_simchiAB(ZERO);
    initV_BB.set_chiS(6.107227423471418e-04);
    initV_BB.set_chiT(3.244360647922483e-04);
    initV_BB.set_chiSPrim(3.244360647922483e-04);
    initV_BB.set_chiTPrim(1.598368940936969e-02);
    initV_BB.set_DeltaS(4.672072595867999e-02);
    initV_BB.set_DeltaT(2.175818210965474e-03);
    initV_BB.set_simDeltaB(ZERO);
    initV_BB.set_mu(5.391034757927860e+00);
    initV_BB.set_mAF(6.509518143510621e-01);
    initV_BB.set_d(8.233828080112927e-02);

//    sdelta                    barchiAB             simchiAB                 chiS                 chiT                     chiSPrim                chiTPrim               DeltaS                  DeltaT                simDeltaB               mu                     mAF                    d                      E                      gt               gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                      gBDelta               gm                 lbarchiAB                 lsimchiAB             lDeltaS              lchiS              lchiT                  lchiSPrim              lchiTPrim              lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
//    0                   2.026312528393326e-01  4.227990941756281e-25 6.891768026409183e-04  8.628367873787484e-14 -3.783852286645283e-04 -4.073204472344117e-14 -3.552596135747061e-18  6.609526106245610e-17  2.857903213852130e-13 5.563039518274006e+00 -4.228220586641046e-13  7.975516757789704e-02   -0.176758926698306     0.057751449458077     0.057751449458074     0.057751449458080     0.057751449458074     0.057751449458080     3.883123799303418     0.057751449458077     1.970564335235827     0.254462009553847     0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000    -0.000000000000000    -0.000000000000092   -11.127243656127373    -0.000000000000555
    // 0
    Var initV_CC(N_VAR);
    initV_CC.set_barchiAB(2.026312528393326e-01);
    initV_CC.set_simchiAB(ZERO);
    initV_CC.set_chiS(6.891768026409183e-04);
    initV_CC.set_chiT(ZERO);
    initV_CC.set_chiSPrim(3.783852286645283e-04);
    initV_CC.set_chiTPrim(ZERO);
    initV_CC.set_DeltaS(ZERO);
    initV_CC.set_DeltaT(ZERO);
    initV_CC.set_simDeltaB(ZERO);
    initV_CC.set_mu(5.563039518274006e+00);
    initV_CC.set_mAF(ZERO);
    initV_CC.set_d(7.975516757789704e-02);

    // DS+DT+m?
    Var initV_DD(N_VAR);
    initV_DD.set_barchiAB(1.715797881677855e-01);
    initV_DD.set_simchiAB(ZERO);
    initV_DD.set_chiS(1.632058553614935e-01);
    initV_DD.set_chiT(6.710615038101641e-05);
    initV_DD.set_chiSPrim(2.477049655097007e-04);
    initV_DD.set_chiTPrim(ZERO);
    initV_DD.set_DeltaS(1.670753025582039e-01);
    initV_DD.set_DeltaT(ZERO);
    initV_DD.set_simDeltaB(ZERO);
    initV_DD.set_mu(3.460201913896650e+00);
    initV_DD.set_mAF(5.365718830046698e-02);
    initV_DD.set_d(2.168172760566418e-02);

    /////////////////////////

    // ZERO - t'=0, t''=0 - do porownania
    // PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, t'=0.000000, t''=0.000000, U=12.000000, beta=500.000000, L=512
    //                    barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    //DS       1.736268397554770e-01  4.560658295008801e-23 5.623331999026073e-04 -1.291482207909358e-13  4.694260638473811e-04 -5.767262310259065e-13  1.647400197495271e-01  6.472704823168465e-13 -1.427599752058793e-18 3.503826560214176e+00  1.075962722721105e-11 3.133850674127103e-02   -0.243326240325883     0.020652813272617     0.020652813272786     0.020652813272447     0.020652813272786     0.020652813272447     3.952555135593992     0.020652813272617     1.988103401635336     0.192220227560733     0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000     0.000000000000000     0.162786002774712    -0.000000000000213     0.000000000000000    -7.014019036553054     0.000000000008679
    //DS+m     1.733802391311728e-01  1.180032036470138e-27 5.917142216861124e-04 -7.201427264640855e-04  3.912035940240002e-04 -7.751836478760639e-03 -1.643694748352157e-01 -8.105564191823008e-05 -7.229039097759481e-21 3.489032051212920e+00  6.592867525937558e-02 3.126274304614424e-02   -0.243337347753924     0.020595761917676     0.021659594031227     0.019584180957318     0.021659594031227     0.019584180957318     3.918491836727725     0.020621887494272     1.979518081940077     0.190443024839027     0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000     0.000000000000000    -0.161020111337254     0.000026467989265     0.000000000000000    -6.984734806255074     0.052223552145509

    Var initV_Z1(N_VAR);
    initV_Z1.set_barchiAB(1.736268397554770e-01);
    initV_Z1.set_simchiAB(ZERO);
    initV_Z1.set_chiS(5.623331999026073e-04);
    initV_Z1.set_chiT(ZERO);
    initV_Z1.set_chiSPrim(4.694260638473811e-04);
    initV_Z1.set_chiTPrim(ZERO);
    initV_Z1.set_DeltaS(1.647400197495271e-01);
    initV_Z1.set_DeltaT(ZERO);
    initV_Z1.set_simDeltaB(ZERO);
    initV_Z1.set_mu(3.503826560214176e+00);
    initV_Z1.set_mAF(ZERO);
    initV_Z1.set_d(3.133850674127103e-02);

    Var initV_Z2(N_VAR);
    initV_Z2.set_barchiAB(1.733802391311728e-01);
    initV_Z2.set_simchiAB(ZERO);
    initV_Z2.set_chiS(5.917142216861124e-04);
    initV_Z2.set_chiT(7.201427264640855e-04);
    initV_Z2.set_chiSPrim(3.912035940240002e-04);
    initV_Z2.set_chiTPrim(7.751836478760639e-03);
    initV_Z2.set_DeltaS(1.643694748352157e-01);
    initV_Z2.set_DeltaT(8.105564191823008e-05);
    initV_Z2.set_simDeltaB(ZERO);
    initV_Z2.set_mu(3.489032051212920e+00);
    initV_Z2.set_mAF(6.592867525937558e-02);
    initV_Z2.set_d(3.126274304614424e-02);

    //50.000000  1.734958478949029e-01 1.479410839312346e-28 5.694947588986846e-04 -4.791833386978826e-04 4.497952965055218e-04 -5.043349909200833e-03 1.645495503298624e-01 5.135141306044836e-05 7.464711704205798e-29 3.496602270064651e+00 4.377315460371044e-02 3.129394352734314e-02   -0.243364556483997     0.020621065069534     0.021321554727265     0.019943589013150     0.021321554727265     0.019943589013150     3.937495645747600     0.020632571870208     1.984312386129663     0.191405851479901     0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000     0.000000000000000     0.161978284483390    -0.000016849663777    -0.000000000000000    -6.999640750939276     0.035049408874484
    Var initV_Z2_B50(N_VAR);
    initV_Z2_B50.set_barchiAB(1.734958478949029e-01);
    initV_Z2_B50.set_simchiAB(ZERO);
    initV_Z2_B50.set_chiS(5.694947588986846e-04);
    initV_Z2_B50.set_chiT(-4.791833386978826e-04);
    initV_Z2_B50.set_chiSPrim(4.497952965055218e-04);
    initV_Z2_B50.set_chiTPrim(-5.043349909200833e-03);
    initV_Z2_B50.set_DeltaS(1.645495503298624e-01);
    initV_Z2_B50.set_DeltaT(5.135141306044836e-05);
    initV_Z2_B50.set_simDeltaB(ZERO);
    initV_Z2_B50.set_mu(3.496602270064651e+00);
    initV_Z2_B50.set_mAF(4.377315460371044e-02);
    initV_Z2_B50.set_d(3.129394352734314e-02);

    // STANDARD - t'=0.25, t''=0
    //PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, t'=0.250000, t''=0.000000, U=12.000000, beta=500.000000, L=512
    //                      barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    // DS         1.736126774702611e-01 -4.374535151143151e-23 -7.105665884589442e-04  3.253559670497917e-13  2.370085064596791e-03  4.210151247296573e-12  1.647364994545274e-01  2.435357064581637e-13  3.054348850582745e-21 3.504606218225224e+00 -3.102341566783024e-11  3.143224122455469e-02   -0.243329264265046     0.020707478060625     0.020707478060137     0.020707478061114     0.020707478060137     0.020707478061114     3.952461554088670     0.020707478060626     1.988079866124264     0.192256836311526    -0.000000000000000    -0.005176869515156     0.000000000000000    -0.000000000000000     0.000000000000000     0.162778670162292    -0.000000000000080    -0.000000000000000    -7.027222728294262    -0.000000000025031
    // DS+DT+m    1.733770448794482e-01  8.150135383561053e-20 -6.390883282464091e-04  7.001707899582296e-04  2.172872463005137e-03  7.516737988623041e-03  1.643909148624374e-01 -7.833765657089206e-05  1.379873610715300e-19 3.489514461003490e+00 -6.405636335579069e-02  3.134804099861067e-02   -0.243338270563175     0.020646155998253     0.019662556713935     0.021678958830522     0.019662556713935     0.021678958830522     3.920307359911151     0.020670757772228     1.979976605899966     0.190568982268390     0.000000000000000    -0.005167689443057     0.000000000000000    -0.000000000000000     0.000000000000000     0.161115728359435     0.000025592307634    -0.000000000000000    -6.997883807312957    -0.050805948849911

    Var initV_S1(N_VAR);
    initV_S1.set_barchiAB(1.736126774702611e-01);
    initV_S1.set_simchiAB(ZERO);
    initV_S1.set_chiS(7.105665884589442e-04);
    initV_S1.set_chiT(ZERO);
    initV_S1.set_chiSPrim(2.370085064596791e-03);
    initV_S1.set_chiTPrim(ZERO);
    initV_S1.set_DeltaS(1.647364994545274e-01);
    initV_S1.set_DeltaT(ZERO);
    initV_S1.set_simDeltaB(ZERO);
    initV_S1.set_mu(3.504606218225224e+00);
    initV_S1.set_mAF(ZERO);
    initV_S1.set_d(3.143224122455469e-02);

    Var initV_S2(N_VAR);
    initV_S2.set_barchiAB(1.733770448794482e-01);
    initV_S2.set_simchiAB(ZERO);
    initV_S2.set_chiS(6.390883282464091e-04);
    initV_S2.set_chiT(7.001707899582296e-04);
    initV_S2.set_chiSPrim(2.172872463005137e-03);
    initV_S2.set_chiTPrim(7.516737988623041e-03);
    initV_S2.set_DeltaS(1.643909148624374e-01);
    initV_S2.set_DeltaT(7.833765657089206e-05);
    initV_S2.set_simDeltaB(ZERO);
    initV_S2.set_mu(3.489514461003490e+00);
    initV_S2.set_mAF(6.405636335579069e-02);
    initV_S2.set_d(3.134804099861067e-02);

    //70.000000  1.734344058817622e-01 -9.298139786014975e-30 -6.753322116406877e-04 -5.982764527920221e-04 2.271421522401173e-03 -6.357936976129183e-03 1.644751883285191e-01 6.165035665419414e-05 1.998431873648260e-30 3.493645811028719e+00 5.467090050101556e-02 3.136975507312282e-02   -0.243346486372792     0.020662065650287     0.021540452775789     0.019819497824885     0.021540452775789     0.019819497824885     3.929001648373242     0.020679975300337     1.982170943277406     0.191018082298807    -0.000000000000000    -0.005169993825084     0.000000000000000    -0.000000000000000     0.000000000000000     0.161555821514813    -0.000020185362743    -0.000000000000000    -7.005728731298370     0.043573209224471
    //50.000000  1.734992180479828e-01 1.611905756637644e-29 -6.971858513513328e-04 -4.369703136638173e-04 2.330203559741241e-03 -4.581294535319059e-03 1.645742112515331e-01 4.199587639734937e-05 -1.090306273336562e-29 3.498260239869022e+00 3.987210358502264e-02 3.139059398176804e-02   -0.243367140533457     0.020678367933159     0.021315445367397     0.020060331511212     0.021315445367397     0.020060331511212     3.939965615369503     0.020687888439305     1.984934662745729     0.191573606283796     0.000000000000000    -0.005171972109826     0.000000000000000    -0.000000000000000     0.000000000000000     0.162104183376899    -0.000013788525749     0.000000000000000    -7.014697192662999     0.031978928404410
    Var initV_S2_B50(N_VAR);
    initV_S2_B50.set_barchiAB(1.734992180479828e-01);
    initV_S2_B50.set_simchiAB(ZERO);
    initV_S2_B50.set_chiS(-6.971858513513328e-04);
    initV_S2_B50.set_chiT(-4.369703136638173e-04);
    initV_S2_B50.set_chiSPrim(2.330203559741241e-03);
    initV_S2_B50.set_chiTPrim(-4.581294535319059e-03);
    initV_S2_B50.set_DeltaS(1.645742112515331e-01);
    initV_S2_B50.set_DeltaT(4.199587639734937e-05);
    initV_S2_B50.set_simDeltaB(ZERO);
    initV_S2_B50.set_mu(3.498260239869022e+00);
    initV_S2_B50.set_mAF(3.987210358502264e-02);
    initV_S2_B50.set_d(3.139059398176804e-02);

    // Beta 50, n=0.9999
    //0.999000  1.693632212068261e-01 3.878155224245762e-33 -1.744746434721240e-04 -3.773764794559694e-04 5.209508357091014e-04 -1.385442962588440e-02 1.669608731112841e-01 3.549990794519945e-05 -4.607401607028525e-32 3.374481699234608e+00 1.284268878610062e-01 1.498401845089659e-02   -0.233049174048928     0.004976367975940     0.005514414069568     0.004490819499505     0.005514414069568     0.004490819499505     3.860042880532571     0.005002616784536     1.964699183216751     0.168413692036808     0.000000000000000    -0.001250654196134     0.000000000000000    -0.000000000000000     0.000000000000000     0.161119032395179    -0.000011419263910     0.000000000000000    -6.753399596496808     0.098652359602822
    Var initV_S2_n0999_B50(N_VAR);
    initV_S2_n0999_B50.set_barchiAB(1.693632212068261e-01);
    initV_S2_n0999_B50.set_simchiAB(ZERO);
    initV_S2_n0999_B50.set_chiS(-1.744746434721240e-04);
    initV_S2_n0999_B50.set_chiT(-3.773764794559694e-04);
    initV_S2_n0999_B50.set_chiSPrim(5.209508357091014e-04);
    initV_S2_n0999_B50.set_chiTPrim(-1.385442962588440e-02);
    initV_S2_n0999_B50.set_DeltaS(1.669608731112841e-01);
    initV_S2_n0999_B50.set_DeltaT(3.549990794519945e-05);
    initV_S2_n0999_B50.set_simDeltaB(ZERO);
    initV_S2_n0999_B50.set_mu(3.374481699234608e+00);
    initV_S2_n0999_B50.set_mAF(1.284268878610062e-01);
    initV_S2_n0999_B50.set_d(1.498401845089659e-02);

    //PARAMETERS: n=0.999000, J=0.333333, t=-1.000000, t'=0.250000, t''=0.000000, U=12.000000, beta=50.000000, L=512
    //#n               barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    //20.000000  1.693213376848892e-01 -6.249758098372775e-32 -1.972878336756076e-04 -2.998863639715042e-04 5.717324674911084e-04 -9.747217165171969e-03 1.669386600807075e-01 1.883745233601774e-05 -2.891887508584630e-32 3.382852699779633e+00 1.006375211430579e-01 1.493953521070092e-02   -0.233581969551131     0.004967936808698     0.005385090574426     0.004583097683152     0.005385090574426     0.004583097683152     3.908866968310579     0.004984094128789     1.977085473192947     0.170431582785354    -0.000000000000000    -0.001246023532197     0.000000000000000    -0.000000000000000     0.000000000000000     0.163135253530876    -0.000006136091267     0.000000000000000    -6.769869843092618     0.080467756348649
    //10.000000  1.608501251110715e-01 3.737104622904047e-32 -2.967376654374462e-04 -6.849339846205974e-04 5.106720601875112e-04 -1.300635558412042e-02 1.576484317600355e-01 2.951954452560371e-06 1.216754861431882e-30 3.082301346878736e+00 2.073494484425890e-01 1.340865656687066e-02   -0.239249706479386     0.004525934062384     0.005428828307023     0.003773204451972     0.005428828307023     0.003773204451972     3.667040043442593     0.004601016379497     1.914951707861739     0.151986896506147     0.000000000000000    -0.001150254094874     0.000000000000000    -0.000000000000000     0.000000000000000     0.144525778012494    -0.000000902077932    -0.000000000000000    -6.168113931062456     0.168755879679064

    //20.000000  1.693213376848892e-01 -6.249758098372775e-32 -1.972878336756076e-04 -2.998863639715042e-04 5.717324674911084e-04 -9.747217165171969e-03 1.669386600807075e-01 1.883745233601774e-05 -2.891887508584630e-32 3.382852699779633e+00 1.006375211430579e-01 1.493953521070092e-02   -0.233581969551131     0.004967936808698     0.005385090574426     0.004583097683152     0.005385090574426     0.004583097683152     3.908866968310579     0.004984094128789     1.977085473192947     0.170431582785354    -0.000000000000000    -0.001246023532197     0.000000000000000    -0.000000000000000     0.000000000000000     0.163135253530876    -0.000006136091267     0.000000000000000    -6.769869843092618     0.080467756348649
    Var initV_S2_n0999_B20(N_VAR);
    initV_S2_n0999_B20.set_barchiAB(1.693213376848892e-01);
    initV_S2_n0999_B20.set_simchiAB(ZERO);
    initV_S2_n0999_B20.set_chiS(-1.972878336756076e-04);
    initV_S2_n0999_B20.set_chiT(-2.998863639715042e-04);
    initV_S2_n0999_B20.set_chiSPrim(5.717324674911084e-04);
    initV_S2_n0999_B20.set_chiTPrim(-9.747217165171969e-03);
    initV_S2_n0999_B20.set_DeltaS(1.669386600807075e-01);
    initV_S2_n0999_B20.set_DeltaT(1.883745233601774e-05);
    initV_S2_n0999_B20.set_simDeltaB(ZERO);
    initV_S2_n0999_B20.set_mu(3.382852699779633e+00);
    initV_S2_n0999_B20.set_mAF(1.006375211430579e-01);
    initV_S2_n0999_B20.set_d(1.493953521070092e-02);

    Var initV_S2_n0999_B10(N_VAR);
    initV_S2_n0999_B10.set_barchiAB(1.608501251110715e-01);
    initV_S2_n0999_B10.set_simchiAB(ZERO);
    initV_S2_n0999_B10.set_chiS(-2.967376654374462e-04);
    initV_S2_n0999_B10.set_chiT(-6.849339846205974e-04);
    initV_S2_n0999_B10.set_chiSPrim(5.106720601875112e-04);
    initV_S2_n0999_B10.set_chiTPrim(-1.300635558412042e-02);
    initV_S2_n0999_B10.set_DeltaS(1.576484317600355e-01);
    initV_S2_n0999_B10.set_DeltaT(2.951954452560371e-06);
    initV_S2_n0999_B10.set_simDeltaB(ZERO);
    initV_S2_n0999_B10.set_mu(3.082301346878736e+00);
    initV_S2_n0999_B10.set_mAF(2.073494484425890e-01);
    initV_S2_n0999_B10.set_d(1.340865656687066e-02);

    // YBCO - t'=0.3, t''=-0.15
    //PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, t'=0.300000, t''=-0.150000, U=12.000000, beta=500.000000, L=512
    //                          barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    //    DS          1.736686634323451e-01  3.406350300256585e-25 -2.160711036653750e-03 -6.343063933539720e-13  5.268720209362227e-03 -5.983698992771900e-12 -1.645906271953541e-01  4.609834128863203e-14 -9.412012034081712e-18 3.559768728361104e+00  5.126802197212828e-11  3.220453198209966e-02   -0.243436769260498     0.021162115695339     0.021162115696146     0.021162115694533     0.021162115696146     0.021162115694533     3.951679929384055     0.021162115695339     1.987883278611714     0.192732858607477     0.000000000000000    -0.006348634708602     0.000000000000000     0.003174317354301     0.000000000000000    -0.162602369513154    -0.000000000000015     0.000000000000000    -7.134267880362792     0.000000000041427
    //    DS+DT+m     1.735292382157748e-01  1.156405778373004e-24 -2.070338749307986e-03 -5.355139016103090e-04  5.021970334178204e-03 -5.595960889538600e-03  1.644030688755860e-01  5.358278227042465e-05  4.129014371802264e-18 3.548373793110221e+00  4.834186287402074e-02  3.211668059177052e-02   -0.243438067624266     0.021103132408608     0.021877458794804     0.020356212375132     0.021877458794804     0.020356212375132     3.933363981455154     0.021116835584968     1.983271030760837     0.191741546240428     0.000000000000000    -0.006335050675490     0.000000000000000     0.003167525337745     0.000000000000000     0.161664277388980    -0.000017563382151    -0.000000000000000    -7.112541149877456     0.038682314550747

    Var initV_Y1(N_VAR);
    initV_Y1.set_barchiAB(1.736686634323451e-01);
    initV_Y1.set_simchiAB(ZERO);
    initV_Y1.set_chiS(2.160711036653750e-03);
    initV_Y1.set_chiT(ZERO);
    initV_Y1.set_chiSPrim(5.268720209362227e-03);
    initV_Y1.set_chiTPrim(ZERO);
    initV_Y1.set_DeltaS(1.645906271953541e-01);
    initV_Y1.set_DeltaT(ZERO);
    initV_Y1.set_simDeltaB(ZERO);
    initV_Y1.set_mu(3.559768728361104e+00);
    initV_Y1.set_mAF(ZERO);
    initV_Y1.set_d(3.220453198209966e-02);

    Var initV_Y2(N_VAR);
    initV_Y2.set_barchiAB(1.735292382157748e-01);
    initV_Y2.set_simchiAB(ZERO);
    initV_Y2.set_chiS(2.070338749307986e-03);
    initV_Y2.set_chiT(5.355139016103090e-04);
    initV_Y2.set_chiSPrim(5.021970334178204e-03);
    initV_Y2.set_chiTPrim(5.595960889538600e-03);
    initV_Y2.set_DeltaS(1.644030688755860e-01);
    initV_Y2.set_DeltaT(5.358278227042465e-05);
    initV_Y2.set_simDeltaB(ZERO);
    initV_Y2.set_mu(3.548373793110221e+00);
    initV_Y2.set_mAF(4.834186287402074e-02);
    initV_Y2.set_d(3.211668059177052e-02);

    //PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, t'=0.300000, t''=-0.150000, U=12.000000, beta=1500.000000, L=256
    //beta                   barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    //1500.000000  1.735272227556184e-01 2.548410076024401e-24 -2.068204078979807e-03 -5.393294871565949e-04 5.016177707295835e-03 -5.637588204908010e-03 1.644005210479826e-01 5.415970531516095e-05 1.477311162013282e-22 3.548181194904051e+00 4.868886682601158e-02 3.211507211957156e-02   -0.243438029928446     0.021102081610089     0.021882066254762     0.020349899460795     0.021882066254762     0.020349899460795     3.933101370309981     0.021115982857779     1.983204823085598     0.191727121011641     0.000000000000000    -0.006334794857334     0.000000000000000     0.003167397428667     0.000000000000000     0.161650978653374    -0.000017751300933    -0.000000000000000    -7.112180232685462     0.038954405705136

    // 100.000000  1.735579033489476e-01 -9.420709537667491e-23 -2.104590301050614e-03 -4.738390511759797e-04 5.114314886068155e-03 -4.926146193404226e-03 1.644387576886656e-01 4.304514597796458e-05 1.051856478317989e-24 3.551299821954484e+00 4.273329854607272e-02 3.214094685301862e-02   -0.243442535245249     0.021119002236076     0.021802075136161     0.020457330445010     0.021802075136161     0.020457330445010     3.937351367343184     0.021129702790586     1.984276031035799     0.191958614252125    -0.000000000000000    -0.006338910837176     0.000000000000000     0.003169455418588     0.000000000000000     0.161863291857421    -0.000014123655364    -0.000000000000000    -7.118004855339835     0.034270977371526
    Var initV_Y2_B100(N_VAR);
    initV_Y2_B100.set_barchiAB(1.735579033489476e-01);
    initV_Y2_B100.set_simchiAB(ZERO);
    initV_Y2_B100.set_chiS(-2.104590301050614e-03);
    initV_Y2_B100.set_chiT(-4.738390511759797e-04);
    initV_Y2_B100.set_chiSPrim(5.114314886068155e-03);
    initV_Y2_B100.set_chiTPrim(-4.926146193404226e-03);
    initV_Y2_B100.set_DeltaS(1.644387576886656e-01);
    initV_Y2_B100.set_DeltaT(4.304514597796458e-05);
    initV_Y2_B100.set_simDeltaB(ZERO);
    initV_Y2_B100.set_mu(3.551299821954484e+00);
    initV_Y2_B100.set_mAF(4.273329854607272e-02);
    initV_Y2_B100.set_d(3.214094685301862e-02);

    //  70.000000  1.735869330616300e-01 5.834814080153588e-30 -2.129453047179878e-03 -3.925101783194780e-04 5.180821427544498e-03 -4.056643696318179e-03 1.644779309837663e-01 3.230821462007322e-05 1.344859953717067e-30 3.554015230366038e+00 3.536453725242941e-02 3.216086466736415e-02   -0.243450726228628     0.021132498964418     0.021696252476739     0.020583393973679     0.021696252476739     0.020583393973679     3.941858530234936     0.021139823225209     1.985411425935425     0.192196282170995     0.000000000000000    -0.006341946967563     0.000000000000000     0.003170973483781     0.000000000000000     0.162087183820938    -0.000010612867616    -0.000000000000000    -7.123146626732560     0.028434513056597
    Var initV_Y2_B70(N_VAR);
    initV_Y2_B70.set_barchiAB(1.735869330616300e-01);
    initV_Y2_B70.set_simchiAB(ZERO);
    initV_Y2_B70.set_chiS(-2.129453047179878e-03);
    initV_Y2_B70.set_chiT(-3.925101783194780e-04);
    initV_Y2_B70.set_chiSPrim(5.180821427544498e-03);
    initV_Y2_B70.set_chiTPrim(-4.056643696318179e-03);
    initV_Y2_B70.set_DeltaS(1.644779309837663e-01);
    initV_Y2_B70.set_DeltaT(3.230821462007322e-05);
    initV_Y2_B70.set_simDeltaB(ZERO);
    initV_Y2_B70.set_mu(3.554015230366038e+00);
    initV_Y2_B70.set_mAF(3.536453725242941e-02);
    initV_Y2_B70.set_d(3.216086466736415e-02);

    //  50.000000  1.736406558015982e-01 2.221400212387857e-29 -2.160329889390379e-03 -9.774708427137353e-05 5.262190079961294e-03 -9.978886682728231e-04 1.645552303951194e-01 6.807123489455918e-06 -2.691089762170437e-28 3.558697884031622e+00 8.791357373103452e-03 3.219052155573833e-02   -0.243473608228775     0.021153569466321     0.021292345826041     0.021015697604310     0.021292345826041     0.021015697604310     3.951083501630452     0.021154021715175     1.987733257162654     0.192670752053818     0.000000000000000    -0.006346206514553     0.000000000000000     0.003173103257276     0.000000000000000     0.162542863980289    -0.000002241292776     0.000000000000000    -7.132137864989308     0.007106533099985
    Var initV_Y2_B50(N_VAR);
    initV_Y2_B50.set_barchiAB(1.736406558015982e-01);
    initV_Y2_B50.set_simchiAB(ZERO);
    initV_Y2_B50.set_chiS(-2.160329889390379e-03);
    initV_Y2_B50.set_chiT(-9.774708427137353e-05);
    initV_Y2_B50.set_chiSPrim(5.262190079961294e-03);
    initV_Y2_B50.set_chiTPrim(-9.978886682728231e-04);
    initV_Y2_B50.set_DeltaS(1.645552303951194e-01);
    initV_Y2_B50.set_DeltaT(6.807123489455918e-06);
    initV_Y2_B50.set_simDeltaB(ZERO);
    initV_Y2_B50.set_mu(3.558697884031622e+00);
    initV_Y2_B50.set_mAF(8.791357373103452e-03);
    initV_Y2_B50.set_d(3.219052155573833e-02);

    // LSCO - t'=0.1, t''=-0.05
    //PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, t'=0.100000, t''=-0.050000, U=12.000000, beta=500.000000, L=512
    //                       barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    //  DS         1.736248615012602e-01  8.501643260495426e-29 -3.263418355523209e-04  2.641532621481224e-13  2.035148732008065e-03  3.198527026843648e-12 -1.647317269771538e-01 -6.716169833093619e-14 -2.473376341529117e-21 3.507797852212508e+00 -2.445430362979755e-11 3.141420478033410e-02   -0.243334626482279     0.020696950777094     0.020696950776709     0.020696950777479     0.020696950776709     0.020696950777479     3.952479582726079     0.020696950777094     1.988084400302482     0.192259130811438     0.000000000000000    -0.002069695077709     0.000000000000000     0.001034847538855     0.000000000000000    -0.162774696876102     0.000000000000022     0.000000000000000    -7.024676445223550    -0.000000000019728
    //  DS+DT+m    1.733894345418540e-01 -1.113708472771950e-16 -2.562324861021548e-04  7.015958769991259e-04  1.843298139945970e-03  7.530756979197292e-03  1.643838174677944e-01 -7.853170372498459e-05  1.768637296107604e-15 3.492762112457172e+00 -6.416258250727620e-02 3.133121921590688e-02   -0.243343770509419     0.020636309914489     0.019651125378842     0.021670885441774     0.019651125378842     0.021670885441774     3.920217645718592     0.020661005410308     1.979953950403542     0.190567390132525    -0.000000000000000    -0.002066100541031     0.000000000000000     0.001033050270515     0.000000000000000     0.161105085476958     0.000025655114224    -0.000000000000001    -6.995429666860255    -0.050881404188532

    Var initV_L1(N_VAR);
    initV_L1.set_barchiAB(1.736248615012602e-01);
    initV_L1.set_simchiAB(ZERO);
    initV_L1.set_chiS(3.263418355523209e-04);
    initV_L1.set_chiT(ZERO);
    initV_L1.set_chiSPrim(2.035148732008065e-03);
    initV_L1.set_chiTPrim(ZERO);
    initV_L1.set_DeltaS(1.647317269771538e-01);
    initV_L1.set_DeltaT(ZERO);
    initV_L1.set_simDeltaB(ZERO);
    initV_L1.set_mu(3.507797852212508e+00);
    initV_L1.set_mAF(ZERO);
    initV_L1.set_d(3.141420478033410e-02);

    Var initV_L2(N_VAR);
    initV_L2.set_barchiAB(1.733894345418540e-01);
    initV_L2.set_simchiAB(ZERO);
    initV_L2.set_chiS(2.562324861021548e-04);
    initV_L2.set_chiT(7.015958769991259e-04);
    initV_L2.set_chiSPrim(1.843298139945970e-03);
    initV_L2.set_chiTPrim(7.530756979197292e-03);
    initV_L2.set_DeltaS(1.643838174677944e-01);
    initV_L2.set_DeltaT(7.853170372498459e-05);
    initV_L2.set_simDeltaB(ZERO);
    initV_L2.set_mu(3.492762112457172e+00);
    initV_L2.set_mAF(6.416258250727620e-02);
    initV_L2.set_d(3.133121921590688e-02);

//    example_1D_sdelta_run(par, initV_AA);
//    example_1D_sdelta_run(par, initV_BB);
//    example_1D_sdelta_run(par, initV_CC);

//    example_1D_sdelta_run(par, initV_Z1);
//    example_1D_sdelta_run(par, initV_Z2);
//    example_1D_sdelta_run(par, initV_Z2_B50);

//    example_1D_sdelta_run(par, initV_S1);
//    example_1D_sdelta_run(par, initV_S2);
//    example_1D_sdelta_run(par, initV_S2_B50);
//    example_1D_sdelta_run(par, initV_S2_n0999_B20);
    example_1D_sdelta_run(par, initV_S2_n0999_B10);



//    example_1D_sdelta_run(par, initV_Y1);
//    example_1D_sdelta_run(par, initV_Y2);
//    example_1D_sdelta_run(par, initV_Y2_B100);
//    example_1D_sdelta_run(par, initV_Y2_B70);
//    example_1D_sdelta_run(par, initV_Y2_B50);

//    example_1D_sdelta_run(par, initV_L1);
//    example_1D_sdelta_run(par, initV_L2);

//    example_1D_beta_run(par, initV_Y2);
//    example_1D_beta_run(par, initV_S2_n09999_B50);
}

void example_inf() {
    Par par(N_PAR);
    par.set_n(0.996);
    par.set_J(1./3.);
    par.set_t(-1.);
    par.set_tPrim(0.1);
    par.set_tPrim2(-0.05);
    par.set_U(12);
    par.set_beta(500.);
    par.set_L(512);
    par.print_verboten();

	// I'm choosing start point
	Var initV(N_VAR);

	FILE * pFile;
	pFile = fopen(OUTPUT_FILE, "w");

	// INFINITE
	Res res(N_VAR);
	RootGSL root(par, function, res);
	printf("sdelta      "); res.print_label();
	while(1) {
	    initV.set_barchiAB(frand());
	    initV.set_simchiAB(frand());
	    initV.set_chiS(frand());
	    initV.set_chiT(frand());
	    initV.set_chiSPrim(frand());
	    initV.set_chiTPrim(frand());
	    initV.set_DeltaS(frand());
	    initV.set_DeltaT(frand());
	    initV.set_simDeltaB(frand());
	    initV.set_mu(16.*frand()-8.);
	    initV.set_mAF(frand());
	    initV.set_d(frand());

		if(root.Solve(initV) == GSL_SUCCESS) {
			printf("       "); root.print(); root.print(pFile); fflush(pFile); fflush(stdout);
		}
	}
	fclose (pFile);
}

void example_2D_nt_run(Par &par, Var &initV, int s1, int s2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	int LOOPS_t, LOOPS_n;
	Var initV_cp(N_VAR); initV_cp = initV;

	double start_tPrim = 0.;
	double start_tPrim2 = 0.;
	double start_n = 0.996;

	par.set_tPrim(start_tPrim);
	par.set_tPrim2(start_tPrim2);
	par.set_n(start_n);

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#tPrim   tPrim2    sdelta        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	if(s1<0) { LOOPS_t = 0; } else { LOOPS_t = 100; }
	if(s2<0) { LOOPS_n = 124; } else { LOOPS_n = 0; }

	for(int i=0; i<=LOOPS_t; i++) {
		if(s1<0) { par.set_tPrim(start_tPrim - (double)i*0.01); } else { par.set_tPrim(start_tPrim + (double)i*0.01); }
		par.set_tPrim2(-par.tPrim()/2.);
		par.set_n(start_n);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS_n; j++) {
			if(s2<0) { par.set_n(start_n - (double)j*0.004); } else { par.set_n(start_n + (double)j*0.004); }
			printf("%f %f %f", par.tPrim(), par.tPrim2(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
	}
	printf("\n#END\n");
	fflush(stdout);
	initV = initV_cp;
}

void example_2D_nbeta_run(Par &par, Var &initV, int s2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	//int LOOPS_n;
	Var initV_cp(N_VAR); initV_cp = initV;

	double start_beta= 1500.;
	double start_n = 0.999;

	par.set_beta(start_beta);
	par.set_n(start_n);

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#beta      sdelta        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	//int MNOZNIK = 1;
	//if(s2<0) { LOOPS_n = 250*MNOZNIK; } else { LOOPS_n = 4*MNOZNIK; }

	while(par.beta() > 300) {
		par.set_n(start_n);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=17; j++) {
			par.set_n(0.999 - (double)j*0.0005);
			printf("%f %f ", par.beta(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}
		for(int j=0; j<=98; j++) {
			par.set_n(0.99 - (double)j*0.005);
			printf("%f %f ", par.beta(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
		par.set_beta(par.beta() - 200.);
	}

	while(par.beta() > 100) {
		par.set_n(start_n);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=17; j++) {
			par.set_n(0.999 - (double)j*0.0005);
			printf("%f %f ", par.beta(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}
		for(int j=0; j<=98; j++) {
			par.set_n(0.99 - (double)j*0.005);
			printf("%f %f ", par.beta(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
		par.set_beta(par.beta() - 50.);
	}

	while(par.beta() > 50) {
		par.set_n(start_n);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=17; j++) {
			par.set_n(0.999 - (double)j*0.0005);
			printf("%f %f ", par.beta(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}
		for(int j=0; j<=98; j++) {
			par.set_n(0.99 - (double)j*0.005);
			printf("%f %f ", par.beta(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
		par.set_beta(par.beta() - 10.);
	}

	while(par.beta() > 0) {
		par.set_n(start_n);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=17; j++) {
			par.set_n(0.999 - (double)j*0.0005);
			printf("%f %f ", par.beta(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}
		for(int j=0; j<=98; j++) {
			par.set_n(0.99 - (double)j*0.005);
			printf("%f %f ", par.beta(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
		par.set_beta(par.beta() - 1.);
	}
	printf("\n#END\n");
	fflush(stdout);
	initV = initV_cp;
}

void example_2D_betaN_run(Par &par, Var &initV, int s2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	//int LOOPS_n;
	Var initV_cp(N_VAR); initV_cp = initV;

	double start_beta= 1500.;
	double start_n = 0.999;

	par.set_beta(start_beta);
	par.set_n(start_n);

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      beta        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	//int MNOZNIK = 1;
	//if(s2<0) { LOOPS_n = 250*MNOZNIK; } else { LOOPS_n = 4*MNOZNIK; }

	while(par.n() > 0.99) {
		par.set_beta(start_beta);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { par.set_n(par.n() - 0.005); continue; }
		old.set(res.get_V());
		initV = old;

		while(par.beta() > 300) {
			printf("%f %f ", par.n(), par.beta());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); break; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
			par.set_beta(par.beta() - 200.);
		}
		while(par.beta() <= 300 && par.beta() > 100) {
			printf("%f %f ", par.n(), par.beta());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); break; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
			par.set_beta(par.beta() - 50.);
		}
		while(par.beta() <= 100 && par.beta() > 50) {
			printf("%f %f ", par.n(), par.beta());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); break; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
			par.set_beta(par.beta() - 10.);
		}
		while(par.beta() <=50 && par.beta() > 0) {
			printf("%f %f ", par.n(), par.beta());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); break; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
			par.set_beta(par.beta() - 1.);
		}

		printf("\n");
		fflush(stdout);
		par.set_n(par.n() - 0.0005);
	}

	while(par.n() > 0.5) {
		par.set_beta(start_beta);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		while(par.beta() > 300) {
			printf("%f %f ", par.n(), par.beta());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
			par.set_beta(par.beta() - 200.);
		}
		while(par.beta() > 100) {
			printf("%f %f ", par.n(), par.beta());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
			par.set_beta(par.beta() - 50.);
		}
		while(par.beta() > 50) {
			printf("%f %f ", par.n(), par.beta());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
			par.set_beta(par.beta() - 10.);
		}
		while(par.beta() > 0) {
			printf("%f %f ", par.n(), par.beta());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
			par.set_beta(par.beta() - 1.);
		}

		printf("\n");
		fflush(stdout);
		par.set_n(par.n() - 0.005);
	}



	printf("\n#END\n");
	fflush(stdout);
	initV = initV_cp;
}

void example_2D() {
	Par par(N_PAR);
	par.set_n(0.999);
	par.set_J(1./3.);
	par.set_t(-1.);
	par.set_tPrim(0.);
	par.set_tPrim2(-0.);
	par.set_U(12);
	par.set_beta(1500);
	par.set_L(512);
	printf("#"); par.print_verboten(); fflush(stdout);

    double ZERO = 1e-12;

    // STANDARD - t'=0.25, t''=0
    //PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, t'=0.250000, t''=0.000000, U=12.000000, beta=500.000000, L=512
    //                      barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    // DS         1.736126774702611e-01 -4.374535151143151e-23 -7.105665884589442e-04  3.253559670497917e-13  2.370085064596791e-03  4.210151247296573e-12  1.647364994545274e-01  2.435357064581637e-13  3.054348850582745e-21 3.504606218225224e+00 -3.102341566783024e-11  3.143224122455469e-02   -0.243329264265046     0.020707478060625     0.020707478060137     0.020707478061114     0.020707478060137     0.020707478061114     3.952461554088670     0.020707478060626     1.988079866124264     0.192256836311526    -0.000000000000000    -0.005176869515156     0.000000000000000    -0.000000000000000     0.000000000000000     0.162778670162292    -0.000000000000080    -0.000000000000000    -7.027222728294262    -0.000000000025031
    // DS+DT+m    1.733770448794482e-01  8.150135383561053e-20 -6.390883282464091e-04  7.001707899582296e-04  2.172872463005137e-03  7.516737988623041e-03  1.643909148624374e-01 -7.833765657089206e-05  1.379873610715300e-19 3.489514461003490e+00 -6.405636335579069e-02  3.134804099861067e-02   -0.243338270563175     0.020646155998253     0.019662556713935     0.021678958830522     0.019662556713935     0.021678958830522     3.920307359911151     0.020670757772228     1.979976605899966     0.190568982268390     0.000000000000000    -0.005167689443057     0.000000000000000    -0.000000000000000     0.000000000000000     0.161115728359435     0.000025592307634    -0.000000000000000    -6.997883807312957    -0.050805948849911

    Var initV_S2(N_VAR);
    initV_S2.set_barchiAB(1.733770448794482e-01);
    initV_S2.set_simchiAB(ZERO);
    initV_S2.set_chiS(6.390883282464091e-04);
    initV_S2.set_chiT(7.001707899582296e-04);
    initV_S2.set_chiSPrim(2.172872463005137e-03);
    initV_S2.set_chiTPrim(7.516737988623041e-03);
    initV_S2.set_DeltaS(1.643909148624374e-01);
    initV_S2.set_DeltaT(7.833765657089206e-05);
    initV_S2.set_simDeltaB(ZERO);
    initV_S2.set_mu(3.489514461003490e+00);
    initV_S2.set_mAF(6.405636335579069e-02);
    initV_S2.set_d(3.134804099861067e-02);


    // STANDARD VALUES //////////////////////////////////////////////////////////////////////////////////////
    // n=0.996, beta=1500, U=12

    Var initV_Z2_U12_n996_B1500(N_VAR);
    initV_Z2_U12_n996_B1500.set_barchiAB(1.733802391311728e-01);
    initV_Z2_U12_n996_B1500.set_simchiAB(ZERO);
    initV_Z2_U12_n996_B1500.set_chiS(5.917142216861124e-04);
    initV_Z2_U12_n996_B1500.set_chiT(7.201427264640855e-04);
    initV_Z2_U12_n996_B1500.set_chiSPrim(3.912035940240002e-04);
    initV_Z2_U12_n996_B1500.set_chiTPrim(7.751836478760639e-03);
    initV_Z2_U12_n996_B1500.set_DeltaS(1.643694748352157e-01);
    initV_Z2_U12_n996_B1500.set_DeltaT(8.105564191823008e-05);
    initV_Z2_U12_n996_B1500.set_simDeltaB(ZERO);
    initV_Z2_U12_n996_B1500.set_mu(3.489032051212920e+00);
    initV_Z2_U12_n996_B1500.set_mAF(6.592867525937558e-02);
    initV_Z2_U12_n996_B1500.set_d(3.126274304614424e-02);

    Var initV_S2_U12_n996_B1500(N_VAR);
    initV_S2_U12_n996_B1500.set_barchiAB(1.733770448794482e-01);
    initV_S2_U12_n996_B1500.set_simchiAB(ZERO);
    initV_S2_U12_n996_B1500.set_chiS(6.390883282464091e-04);
    initV_S2_U12_n996_B1500.set_chiT(7.001707899582296e-04);
    initV_S2_U12_n996_B1500.set_chiSPrim(2.172872463005137e-03);
    initV_S2_U12_n996_B1500.set_chiTPrim(7.516737988623041e-03);
    initV_S2_U12_n996_B1500.set_DeltaS(1.643909148624374e-01);
    initV_S2_U12_n996_B1500.set_DeltaT(7.833765657089206e-05);
    initV_S2_U12_n996_B1500.set_simDeltaB(ZERO);
    initV_S2_U12_n996_B1500.set_mu(3.489514461003490e+00);
    initV_S2_U12_n996_B1500.set_mAF(6.405636335579069e-02);
    initV_S2_U12_n996_B1500.set_d(3.134804099861067e-02);

    // new: "D" as J=0.333333, t=-1.000000, t'=0.250000, t''=-0.125000
    //                                   barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    //0.250000 -0.125000 0.996000 1.734616462124874e-01 -1.568173616771527e-24 -1.579897120089935e-03 -6.252730468250607e-04 4.137495893068586e-03 -6.612648148888489e-03 1.643804904644774e-01 6.671659020285137e-05 -6.568654607784542e-23 3.526418119682210e+00 5.673771466551041e-02 3.181696699693445e-02   -0.243403555649339     0.020923506507500     0.021835382892358     0.020049711366528     0.021835382892358     0.020049711366528     3.926751692453118     0.020942547129443     1.981603313595615     0.191208709717647    -0.000000000000000    -0.005235636782361     0.000000000000000     0.002617818391180     0.000000000000000     0.161370342284415    -0.000021831623624     0.000000000000000    -7.067406856163103     0.045207547172152
    Var initV_D2_U12_n996_B1500(N_VAR);
    initV_D2_U12_n996_B1500.set_barchiAB(1.734616462124874e-01);
    initV_D2_U12_n996_B1500.set_simchiAB(ZERO);
    initV_D2_U12_n996_B1500.set_chiS(-1.579897120089935e-03);
    initV_D2_U12_n996_B1500.set_chiT(-6.252730468250607e-04);
    initV_D2_U12_n996_B1500.set_chiSPrim(4.137495893068586e-03);
    initV_D2_U12_n996_B1500.set_chiTPrim(-6.612648148888489e-03);
    initV_D2_U12_n996_B1500.set_DeltaS(1.643804904644774e-01);
    initV_D2_U12_n996_B1500.set_DeltaT(6.671659020285137e-05);
    initV_D2_U12_n996_B1500.set_simDeltaB(ZERO);
    initV_D2_U12_n996_B1500.set_mu(3.526418119682210e+00);
    initV_D2_U12_n996_B1500.set_mAF(5.673771466551041e-02);
    initV_D2_U12_n996_B1500.set_d(3.181696699693445e-02);

	//#PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, t'=0.000000, t''=0.000000, U=100.000000, beta=1500.000000, L=512
    //#tPrim   tPrim2    sdelta            barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    //0.000000 -0.000000 0.996000 1.709720669562160e-01 -5.208062358674675e-23 5.583125730143783e-04 2.927289961329167e-04 4.206040167951171e-04 7.457365795560003e-03 1.671398226458470e-01 -7.368652384135338e-05 -3.533532341595976e-23 2.324094382028890e+00 -6.302003833674190e-02 1.813654870955317e-03   -0.238834044837805     0.008408609473603     0.007467826148130     0.009467910992715     0.007467826148130     0.009467910992715     3.936805716526834     0.008467868570422     1.984138532594645     0.176679562113513    -0.000000000000000    -0.000000000000000     0.000000000000000     0.000000000000000     0.000000000000000     0.164499252312863     0.000024174127357     0.000000000000000    -4.654574864956196    -0.048501324131362
    //0.250000 -0.125000 0.996000 1.709730305254404e-01 1.455069777295424e-23 -3.402788372735323e-04 2.847891528188946e-04 1.935697105432674e-03 7.250412280695190e-03 1.671574591988954e-01 -7.145525770563202e-05 -1.795403223364732e-22 2.328713128616330e+00 -6.140775049304689e-02 1.816970877006969e-03   -0.238842155294413     0.008410993193204     0.007492795793787     0.009441710203124     0.007492795793787     0.009441710203124     3.938379326518896     0.008467252998455     1.984535040385756     0.176750155396624     0.000000000000000    -0.002116813249614     0.000000000000000     0.001058406624807     0.000000000000000     0.164582370395589     0.000023451492477     0.000000000000000    -4.666720624422558    -0.047314467179076
    Var initV_D2_U100_n996_B1500(N_VAR);
    initV_D2_U100_n996_B1500.set_barchiAB(1.709730305254404e-01);
    initV_D2_U100_n996_B1500.set_simchiAB(ZERO);
    initV_D2_U100_n996_B1500.set_chiS(-3.402788372735323e-04);
    initV_D2_U100_n996_B1500.set_chiT(2.847891528188946e-04);
    initV_D2_U100_n996_B1500.set_chiSPrim(1.935697105432674e-03);
    initV_D2_U100_n996_B1500.set_chiTPrim(7.250412280695190e-03);
    initV_D2_U100_n996_B1500.set_DeltaS(1.671574591988954e-01);
    initV_D2_U100_n996_B1500.set_DeltaT(-7.145525770563202e-05);
    initV_D2_U100_n996_B1500.set_simDeltaB(ZERO);
    initV_D2_U100_n996_B1500.set_mu(2.328713128616330e+00);
    initV_D2_U100_n996_B1500.set_mAF(-6.140775049304689e-02);
    initV_D2_U100_n996_B1500.set_d(1.816970877006969e-03);

    // VALUES FOR n=0.999
	//#beta      sdelta            barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    //#PARAMETERS: n=_______, J=0.333333, t=-1.000000, t'=0.000000, t''=0.000000, U=12.000000, beta=1500.000000, L=256
    //1500.000000 0.999000  1.693217473268961e-01 7.924341876053711e-39 1.442308124273939e-04 -3.865254347226006e-04 8.791624100731755e-05 -1.434708493837335e-02 1.669293228458277e-01 3.915289953743466e-05 -5.440517417514205e-36 3.371567247648004e+00 1.317210700617173e-01 1.496803286385485e-02   -0.233050319795889     0.004971261856308     0.005523841203971     0.004473959973037     0.005523841203971     0.004473959973037     3.853549916838822     0.004998900588504     1.963046081180679     0.168093713188949     0.000000000000000    -0.000000000000000     0.000000000000000    -0.000000000000000     0.000000000000000     0.160817619542625    -0.000012573137730     0.000000000000000    -6.744776026624029     0.100744811264657
    //#PARAMETERS: n=_______, J=0.333333, t=-1.000000, t'=0.250000, t''=0.000000, U=12.000000, beta=1500.000000, L=256
    //1500.000000 0.999000  1.693209725050858e-01 4.123560470044455e-39 -1.667374994561771e-04 -3.865208109545345e-04 4.989562599005346e-04 -1.434284630582273e-02 1.669300254048142e-01 3.913925867914187e-05 2.555056580402850e-36 3.371759590640090e+00 1.316929293145910e-01 1.497875572281256e-02   -0.233050431833151     0.004974328326067     0.005526765401681     0.004477111021935     0.005526765401681     0.004477111021935     3.853601116359647     0.005001938211808     1.963059121972552     0.168098200493242     0.000000000000000    -0.001250484552952     0.000000000000000    -0.000000000000000     0.000000000000000     0.160820433063484    -0.000012568924245    -0.000000000000000    -6.748046198759600     0.100728491298010

    Var initV_Z2_U12_n999_B1500(N_VAR);
    initV_Z2_U12_n999_B1500.set_barchiAB(1.693217473268961e-01);
    initV_Z2_U12_n999_B1500.set_simchiAB(ZERO);
    initV_Z2_U12_n999_B1500.set_chiS(1.442308124273939e-04);
    initV_Z2_U12_n999_B1500.set_chiT(-3.865254347226006e-04);
    initV_Z2_U12_n999_B1500.set_chiSPrim(8.791624100731755e-05);
    initV_Z2_U12_n999_B1500.set_chiTPrim(-1.434708493837335e-02);
    initV_Z2_U12_n999_B1500.set_DeltaS(1.669293228458277e-01);
    initV_Z2_U12_n999_B1500.set_DeltaT(3.915289953743466e-05);
    initV_Z2_U12_n999_B1500.set_simDeltaB(ZERO);
    initV_Z2_U12_n999_B1500.set_mu(3.371567247648004e+00);
    initV_Z2_U12_n999_B1500.set_mAF(1.317210700617173e-01);
    initV_Z2_U12_n999_B1500.set_d(1.496803286385485e-02);

    Var initV_S2_U12_n999_B1500(N_VAR);
    initV_S2_U12_n999_B1500.set_barchiAB(1.693209725050858e-01);
    initV_S2_U12_n999_B1500.set_simchiAB(ZERO);
    initV_S2_U12_n999_B1500.set_chiS(-1.667374994561771e-04);
    initV_S2_U12_n999_B1500.set_chiT(-3.865208109545345e-04);
    initV_S2_U12_n999_B1500.set_chiSPrim(4.989562599005346e-04);
    initV_S2_U12_n999_B1500.set_chiTPrim(-1.434284630582273e-02);
    initV_S2_U12_n999_B1500.set_DeltaS(1.669300254048142e-01);
    initV_S2_U12_n999_B1500.set_DeltaT(3.913925867914187e-05);
    initV_S2_U12_n999_B1500.set_simDeltaB(ZERO);
    initV_S2_U12_n999_B1500.set_mu(3.371759590640090e+00);
    initV_S2_U12_n999_B1500.set_mAF(1.316929293145910e-01);
    initV_S2_U12_n999_B1500.set_d(1.497875572281256e-02);

	//#beta      sdelta            barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    //#PARAMETERS: n=_______, J=0.333333, t=-1.000000, t'=0.250000, t''=-0.125000, U=12.000000, beta=1500.000000, L=256
    //1500.000000 0.999000  1.704958772844220e-01 -6.171088531278019e-39 -4.323437435306862e-04 1.477142791670935e-19 1.101517563241439e-03 7.423776040119731e-17 1.682272991251121e-01 -5.881247321528863e-18 1.619737213351462e-37 3.444834324705562e+00 -5.984288244045741e-16 1.521416978082985e-02   -0.232932493225433     0.005054086334453     0.005054086334453     0.005054086334453     0.005054086334453     0.005054086334453     3.988313002024394     0.005054086334454     1.997076113227634     0.175051817375707    -0.000000000000000    -0.001263521583613     0.000000000000000     0.000631760791807     0.000000000000000     0.167735781099033     0.000000000000000    -0.000000000000000    -6.892530590727684    -0.000000000000000
    //1500.000000 0.999000  1.693245461182887e-01 2.106768520669785e-39 -3.736238620659123e-04 -3.868801680576318e-04 9.250397210152400e-04 -1.432984573403581e-02 1.669273918827868e-01 3.908514799277455e-05 9.514228240611516e-39 3.379022813803962e+00 1.316023919500209e-01 1.502374526920111e-02   -0.233053996239784     0.004987205499107     0.005539167323731     0.004490245056106     0.005539167323731     0.004490245056106     3.853760878162308     0.005014706189918     1.963099813601516     0.168121283384920     0.000000000000000    -0.001253676547480     0.000000000000000     0.000626838273740     0.000000000000000     0.160824563082888    -0.000012552067854    -0.000000000000000    -6.761697143457932     0.100674690139747
    //#PARAMETERS: n=_______, J=0.333333, t=-1.000000, t'=0.250000, t''=-0.125000, U=100.000000, beta=1500.000000, L=256
    //1500.000000 0.999000  1.686362278267491e-01 -4.649104688979895e-31 -8.216100226479327e-05 1.638747993298522e-04 4.500273716699871e-04 1.434993893051826e-02 1.676355146110538e-01 -3.865334093159104e-05 -1.198587789363651e-30 2.267519175807586e+00 -1.317354966796765e-01 9.086775374475914e-04   -0.231958465027650     0.002084060311582     0.001623416544443     0.002675411555451     0.001623416544443     0.002675411555451     3.856967558468685     0.002149414049947     1.963916382758870     0.164690175289159    -0.000000000000000    -0.000537353512487     0.000000000000000     0.000268676756243     0.000000000000000     0.161641185375509     0.000012423723500     0.000000000000000    -4.537496283350060    -0.099984478357976

    Var initV_D2_U12_n999_B1500(N_VAR);
    initV_D2_U12_n999_B1500.set_barchiAB(1.693245461182887e-01);
    initV_D2_U12_n999_B1500.set_simchiAB(ZERO);
    initV_D2_U12_n999_B1500.set_chiS(-3.736238620659123e-04);
    initV_D2_U12_n999_B1500.set_chiT(-3.868801680576318e-04);
    initV_D2_U12_n999_B1500.set_chiSPrim(9.250397210152400e-04);
    initV_D2_U12_n999_B1500.set_chiTPrim(-1.432984573403581e-02);
    initV_D2_U12_n999_B1500.set_DeltaS(1.669273918827868e-01);
    initV_D2_U12_n999_B1500.set_DeltaT(3.908514799277455e-05);
    initV_D2_U12_n999_B1500.set_simDeltaB(ZERO);
    initV_D2_U12_n999_B1500.set_mu(3.379022813803962e+00);
    initV_D2_U12_n999_B1500.set_mAF(1.316023919500209e-01);
    initV_D2_U12_n999_B1500.set_d(1.502374526920111e-02);

    Var initV_D2_U100_n999_B1500(N_VAR);
    initV_D2_U100_n999_B1500.set_barchiAB(1.686362278267491e-01);
    initV_D2_U100_n999_B1500.set_simchiAB(ZERO);
    initV_D2_U100_n999_B1500.set_chiS(-8.216100226479327e-05);
    initV_D2_U100_n999_B1500.set_chiT(1.638747993298522e-04);
    initV_D2_U100_n999_B1500.set_chiSPrim(4.500273716699871e-04);
    initV_D2_U100_n999_B1500.set_chiTPrim(1.434993893051826e-02);
    initV_D2_U100_n999_B1500.set_DeltaS(1.676355146110538e-01);
    initV_D2_U100_n999_B1500.set_DeltaT(-3.865334093159104e-05);
    initV_D2_U100_n999_B1500.set_simDeltaB(ZERO);
    initV_D2_U100_n999_B1500.set_mu(2.267519175807586e+00);
    initV_D2_U100_n999_B1500.set_mAF(-1.317354966796765e-01);
    initV_D2_U100_n999_B1500.set_d(9.086775374475914e-04);


    //printf("# CASE EJ SGA_v2+-\n"); example_2D_nt_run(par, initV_S2, 1, -1);
    //printf("# CASE nBeta SGA_v2+\n"); example_2D_nbeta_run(par, initV_D2_U12_n996_B1500, 1);

    //printf("# CASE nBeta SGA_v2- (initV_Z2_U12_n999_B1500)\n"); example_2D_nbeta_run(par, initV_Z2_U12_n999_B1500, -1);
    //printf("# CASE nBeta SGA_v2- (initV_S2_U12_n999_B1500)\n"); example_2D_nbeta_run(par, initV_S2_U12_n999_B1500, -1);
    //printf("# CASE nBeta SGA_v2- (initV_D2_U12_n999_B1500)\n"); example_2D_nbeta_run(par, initV_D2_U12_n999_B1500, -1);
    //printf("# CASE nBeta SGA_v2- (initV_D2_U100_n999_B1500)\n"); example_2D_nbeta_run(par, initV_D2_U100_n999_B1500, -1);


    printf("# CASE BetaN SGA_v2- (initV_Z2_U12_n999_B1500)\n"); example_2D_betaN_run(par, initV_Z2_U12_n999_B1500, -1);
    //printf("# CASE BetaN SGA_v2- (initV_S2_U12_n999_B1500)\n"); example_2D_betaN_run(par, initV_S2_U12_n999_B1500, -1);
    //printf("# CASE BetaN SGA_v2- (initV_D2_U12_n999_B1500)\n"); example_2D_betaN_run(par, initV_D2_U12_n999_B1500, -1);
    //printf("# CASE BetaN SGA_v2- (initV_D2_U100_n999_B1500)\n"); example_2D_betaN_run(par, initV_D2_U100_n999_B1500, -1);
}

void example_1D_v2() {
	printf("#PLOT1D - 2D3t_slow\n");

    Par initPar(N_PAR);
    initPar.set_n(0.996);
    initPar.set_J(1./3.);
    initPar.set_t(-1.);
    initPar.set_tPrim(0.25);
    initPar.set_tPrim2(0.);
    initPar.set_U(12);
    initPar.set_beta(1500.);
    initPar.set_L(512);
    initPar.print_verboten();

    double ZERO = 1e-15;

    // STANDARD - t'=0.25, t''=0
    //PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, t'=0.250000, t''=0.000000, U=12.000000, beta=500.000000, L=512
    //                      barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    // DS         1.736126774702611e-01 -4.374535151143151e-23 -7.105665884589442e-04  3.253559670497917e-13  2.370085064596791e-03  4.210151247296573e-12  1.647364994545274e-01  2.435357064581637e-13  3.054348850582745e-21 3.504606218225224e+00 -3.102341566783024e-11  3.143224122455469e-02   -0.243329264265046     0.020707478060625     0.020707478060137     0.020707478061114     0.020707478060137     0.020707478061114     3.952461554088670     0.020707478060626     1.988079866124264     0.192256836311526    -0.000000000000000    -0.005176869515156     0.000000000000000    -0.000000000000000     0.000000000000000     0.162778670162292    -0.000000000000080    -0.000000000000000    -7.027222728294262    -0.000000000025031
    // DS+DT+m    1.733770448794482e-01  8.150135383561053e-20 -6.390883282464091e-04  7.001707899582296e-04  2.172872463005137e-03  7.516737988623041e-03  1.643909148624374e-01 -7.833765657089206e-05  1.379873610715300e-19 3.489514461003490e+00 -6.405636335579069e-02  3.134804099861067e-02   -0.243338270563175     0.020646155998253     0.019662556713935     0.021678958830522     0.019662556713935     0.021678958830522     3.920307359911151     0.020670757772228     1.979976605899966     0.190568982268390     0.000000000000000    -0.005167689443057     0.000000000000000    -0.000000000000000     0.000000000000000     0.161115728359435     0.000025592307634    -0.000000000000000    -6.997883807312957    -0.050805948849911

    Var initV_S1(N_VAR);
    initV_S1.set_barchiAB(1.736126774702611e-01);
    initV_S1.set_simchiAB(ZERO);
    initV_S1.set_chiS(7.105665884589442e-04);
    initV_S1.set_chiT(ZERO);
    initV_S1.set_chiSPrim(2.370085064596791e-03);
    initV_S1.set_chiTPrim(ZERO);
    initV_S1.set_DeltaS(1.647364994545274e-01);
    initV_S1.set_DeltaT(ZERO);
    initV_S1.set_simDeltaB(ZERO);
    initV_S1.set_mu(3.504606218225224e+00);
    initV_S1.set_mAF(ZERO);
    initV_S1.set_d(3.143224122455469e-02);

    Var initV_S2(N_VAR);
    initV_S2.set_barchiAB(1.733770448794482e-01);
    initV_S2.set_simchiAB(ZERO);
    initV_S2.set_chiS(6.390883282464091e-04);
    initV_S2.set_chiT(7.001707899582296e-04);
    initV_S2.set_chiSPrim(2.172872463005137e-03);
    initV_S2.set_chiTPrim(7.516737988623041e-03);
    initV_S2.set_DeltaS(1.643909148624374e-01);
    initV_S2.set_DeltaT(7.833765657089206e-05);
    initV_S2.set_simDeltaB(ZERO);
    initV_S2.set_mu(3.489514461003490e+00);
    initV_S2.set_mAF(6.405636335579069e-02);
    initV_S2.set_d(3.134804099861067e-02);

    // YBCO - t'=0.3, t''=-0.15
    //PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, t'=0.300000, t''=-0.150000, U=12.000000, beta=500.000000, L=512
    //                          barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    //    DS          1.736686634323451e-01  3.406350300256585e-25 -2.160711036653750e-03 -6.343063933539720e-13  5.268720209362227e-03 -5.983698992771900e-12 -1.645906271953541e-01  4.609834128863203e-14 -9.412012034081712e-18 3.559768728361104e+00  5.126802197212828e-11  3.220453198209966e-02   -0.243436769260498     0.021162115695339     0.021162115696146     0.021162115694533     0.021162115696146     0.021162115694533     3.951679929384055     0.021162115695339     1.987883278611714     0.192732858607477     0.000000000000000    -0.006348634708602     0.000000000000000     0.003174317354301     0.000000000000000    -0.162602369513154    -0.000000000000015     0.000000000000000    -7.134267880362792     0.000000000041427
    //    DS+DT+m     1.735292382157748e-01  1.156405778373004e-24 -2.070338749307986e-03 -5.355139016103090e-04  5.021970334178204e-03 -5.595960889538600e-03  1.644030688755860e-01  5.358278227042465e-05  4.129014371802264e-18 3.548373793110221e+00  4.834186287402074e-02  3.211668059177052e-02   -0.243438067624266     0.021103132408608     0.021877458794804     0.020356212375132     0.021877458794804     0.020356212375132     3.933363981455154     0.021116835584968     1.983271030760837     0.191741546240428     0.000000000000000    -0.006335050675490     0.000000000000000     0.003167525337745     0.000000000000000     0.161664277388980    -0.000017563382151    -0.000000000000000    -7.112541149877456     0.038682314550747

    Var initV_Y1(N_VAR);
    initV_Y1.set_barchiAB(1.736686634323451e-01);
    initV_Y1.set_simchiAB(ZERO);
    initV_Y1.set_chiS(2.160711036653750e-03);
    initV_Y1.set_chiT(ZERO);
    initV_Y1.set_chiSPrim(5.268720209362227e-03);
    initV_Y1.set_chiTPrim(ZERO);
    initV_Y1.set_DeltaS(1.645906271953541e-01);
    initV_Y1.set_DeltaT(ZERO);
    initV_Y1.set_simDeltaB(ZERO);
    initV_Y1.set_mu(3.559768728361104e+00);
    initV_Y1.set_mAF(ZERO);
    initV_Y1.set_d(3.220453198209966e-02);

    Var initV_Y2(N_VAR);
    initV_Y2.set_barchiAB(1.735292382157748e-01);
    initV_Y2.set_simchiAB(ZERO);
    initV_Y2.set_chiS(2.070338749307986e-03);
    initV_Y2.set_chiT(5.355139016103090e-04);
    initV_Y2.set_chiSPrim(5.021970334178204e-03);
    initV_Y2.set_chiTPrim(5.595960889538600e-03);
    initV_Y2.set_DeltaS(1.644030688755860e-01);
    initV_Y2.set_DeltaT(5.358278227042465e-05);
    initV_Y2.set_simDeltaB(ZERO);
    initV_Y2.set_mu(3.548373793110221e+00);
    initV_Y2.set_mAF(4.834186287402074e-02);
    initV_Y2.set_d(3.211668059177052e-02);

    // LSCO - t'=0.1, t''=-0.05
    //PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, t'=0.100000, t''=-0.050000, U=12.000000, beta=500.000000, L=512
    //                       barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    //  DS         1.736248615012602e-01  8.501643260495426e-29 -3.263418355523209e-04  2.641532621481224e-13  2.035148732008065e-03  3.198527026843648e-12 -1.647317269771538e-01 -6.716169833093619e-14 -2.473376341529117e-21 3.507797852212508e+00 -2.445430362979755e-11 3.141420478033410e-02   -0.243334626482279     0.020696950777094     0.020696950776709     0.020696950777479     0.020696950776709     0.020696950777479     3.952479582726079     0.020696950777094     1.988084400302482     0.192259130811438     0.000000000000000    -0.002069695077709     0.000000000000000     0.001034847538855     0.000000000000000    -0.162774696876102     0.000000000000022     0.000000000000000    -7.024676445223550    -0.000000000019728
    //  DS+DT+m    1.733894345418540e-01 -1.113708472771950e-16 -2.562324861021548e-04  7.015958769991259e-04  1.843298139945970e-03  7.530756979197292e-03  1.643838174677944e-01 -7.853170372498459e-05  1.768637296107604e-15 3.492762112457172e+00 -6.416258250727620e-02 3.133121921590688e-02   -0.243343770509419     0.020636309914489     0.019651125378842     0.021670885441774     0.019651125378842     0.021670885441774     3.920217645718592     0.020661005410308     1.979953950403542     0.190567390132525    -0.000000000000000    -0.002066100541031     0.000000000000000     0.001033050270515     0.000000000000000     0.161105085476958     0.000025655114224    -0.000000000000001    -6.995429666860255    -0.050881404188532

    Var initV_L1(N_VAR);
    initV_L1.set_barchiAB(1.736248615012602e-01);
    initV_L1.set_simchiAB(ZERO);
    initV_L1.set_chiS(3.263418355523209e-04);
    initV_L1.set_chiT(ZERO);
    initV_L1.set_chiSPrim(2.035148732008065e-03);
    initV_L1.set_chiTPrim(ZERO);
    initV_L1.set_DeltaS(1.647317269771538e-01);
    initV_L1.set_DeltaT(ZERO);
    initV_L1.set_simDeltaB(ZERO);
    initV_L1.set_mu(3.507797852212508e+00);
    initV_L1.set_mAF(ZERO);
    initV_L1.set_d(3.141420478033410e-02);

    Var initV_L2(N_VAR);
    initV_L2.set_barchiAB(1.733894345418540e-01);
    initV_L2.set_simchiAB(ZERO);
    initV_L2.set_chiS(2.562324861021548e-04);
    initV_L2.set_chiT(7.015958769991259e-04);
    initV_L2.set_chiSPrim(1.843298139945970e-03);
    initV_L2.set_chiTPrim(7.530756979197292e-03);
    initV_L2.set_DeltaS(1.643838174677944e-01);
    initV_L2.set_DeltaT(7.853170372498459e-05);
    initV_L2.set_simDeltaB(ZERO);
    initV_L2.set_mu(3.492762112457172e+00);
    initV_L2.set_mAF(6.416258250727620e-02);
    initV_L2.set_d(3.133121921590688e-02);


	// I'm choosing start point
	Var initV(N_VAR); initV = initV_S2;

	printf("#dochodze do punktu startowego\n");
    automatic_run_J(initPar, initV, 0.33, 0.0, 0.01); printf("\n");
    automatic_run_n(initPar, initV, 0.996, 1.0, 0.001); printf("\n");
    automatic_run_n(initPar, initV, 0.99, 0.6, 0.01); printf("\n");
    //automatic_run_tPrim2(initPar, initV, 0., 0.13, 0.02); printf("\n");

	printf("#Jestem w punkcie startowym\n");

	Par oldPar = initPar;
    Var oldVar = initV;

    automatic_run_tPrim(initPar, initV, 0.25, 0.0, 0.005); initPar=oldPar; initV=oldVar;
    automatic_run_tPrim2(initPar, initV, 0.0, 0.25, 0.005);
    //automatic_run_tPrimANDtPrim2(initPar, initV, 0.25, 0.0, 0.002, 0.5);

    return;
}


void automatic_run_n(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_n, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step+0.5);
	if(LOOPS < 0) { LOOPS = -LOOPS+1; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_n(start + ((double)i)*step);
		printf("%f ", inpar.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	inpar.set_n(stop);
	invar = initV;

	printf("#end automatic_run_n\n\n");
	fflush(stdout);
	return;
}

void automatic_run_tPrim(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_tPrim, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#tPrim       "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step+0.5);
	if(LOOPS < 0) { LOOPS = -LOOPS+1; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_tPrim(start + ((double)i)*step);
		printf("%f ", inpar.tPrim());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_tPrim\n\n");
	fflush(stdout);

	inpar.set_tPrim(stop);
	invar = initV;
	return;
}

void automatic_run_tPrim2(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_tPrim2, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#tPrim2       "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step+0.5);
	if(LOOPS < 0) { LOOPS = -LOOPS+1; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_tPrim2(start + ((double)i)*step);
		printf("%f ", inpar.tPrim2());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_tPrim2\n\n");
	fflush(stdout);

	inpar.set_tPrim2(stop);
	invar = initV;
	return;
}

void automatic_run_tPrimANDtPrim2(Par &inpar, Var &invar, double start, double stop, double step, double ratio) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_tPrimANDtPrim2 | tPrim start %f, stop %f, step %f, ratio %f", start, stop, step, ratio); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#tPrim1  tPrim2        "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step+0.5);
	if(LOOPS < 0) { LOOPS = -LOOPS+1; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_tPrim(start + ((double)i)*step);
		inpar.set_tPrim2(inpar.tPrim()*ratio);
		printf("%f %f ", inpar.tPrim(), inpar.tPrim2());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_tPrimANDtPrim2\n\n");
	fflush(stdout);

	inpar.set_tPrim2(stop);
	invar = initV;
	return;
}

void automatic_run_J(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_J, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#J           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step+0.5);
	if(LOOPS < 0) { LOOPS = -LOOPS+1; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_J(start + ((double)i)*step);
		printf("%f ", inpar.J());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	inpar.set_J(stop);
	invar = initV;

	printf("#end automatic_run_n\n\n");
	fflush(stdout);
	return;
}

} // end of namespace case2D2t_SGA_tPrim0
